\documentclass[twoside,letterpaper,11pt]{article}

%% LaTeX - Article customise

%%% PACKAGES
%\usepackage{booktabs} % for much better looking tables
%\usepackage{array} % for better arrays (eg matrices) in maths
%\usepackage{paralist} % very flexible & customisable lists (eg. enumerate/itemize, etc.)
\usepackage{verbatim} % adds environment for commenting out blocks of text & for better verbatim
%\usepackage{subfigure} % make it possible to include more than one captioned figure/table in a single float
% These packages are all incorporated in the memoir class to one degree or another...

%%% PAGE DIMENSIONS
\usepackage{geometry} % to change the page dimensions
\geometry{margin=0.8in} % for example, change the margins to 2 inches all round
%\geometry{landscape} % set up the page for landscape
% read geometry.pdf for detailed page layout information

%%% HEADERS & FOOTERS
\usepackage{fancyhdr} % This should be set AFTER setting up the page geometry
\pagestyle{fancy} % options: empty , plain , fancy
\renewcommand{\headrulewidth}{0pt} % customise the layout...
\lhead{}\chead{}\rhead{}
\lfoot{}\cfoot{\thepage}\rfoot{}

%%% SECTION TITLE APPEARANCE
%\usepackage{sectsty}
%\allsectionsfont{\sffamily\mdseries\upshape} % (See the fntguide.pdf for font help)
% (This matches ConTeXt defaults)

%%% ToC APPEARANCE
\usepackage[nottoc,notlof,notlot]{tocbibind} % Put the bibliography in the ToC
\usepackage[titles]{tocloft} % Alter the style of the Table of Contents
\renewcommand{\cftsecfont}{\rmfamily\mdseries\upshape}
\renewcommand{\cftsecpagefont}{\rmfamily\mdseries\upshape} % No bold!

%% END Article customise

\usepackage{amsmath}
\usepackage{times}
\renewcommand{\tt}[1]{{\texttt {#1}}}
%\renewcommand{\tt}[1]{\rm #1}
\newcommand{\tmvversion}{0.60}
\newcommand{\myemail}{mike\_jarvis@users.sourceforge.net}
%\newcommand{\website}{http://www.hep.upenn.edu/$\sim$mjarvis/tmv/}
\newcommand{\website}{http://sourceforge.net/projects/tmv-cpp/}

\begin{document}

\title{Template Matrix/Vector Library for C++ \\ User Manual \\ Version \tmvversion}
\author{Mike Jarvis}
\date{April 14, 2007}
\maketitle

%\tableofcontents

\newpage
\section{Overview}

First, this library is provided without any warranty of any kind.  There is no guarantee
that the code will produce accurate results for all inputs.  If someone dies because
my code gave you a wrong result, \underline{do not} blame me.

OK, that was mostly for the crazies out there.  Really, I think this is a pretty good 
product, and I've been using it for my own research extensively.  So at least 
for the routines that I use, they are probably reasonably well debugged.
I also have a test suite which tries to be comprehensive, although
occasionally I find bugs that I hadn't thought to test for in the test suite, so 
there may still be a few lurking in there.  That means the code should definitely be considered
a beta-type release.  Hence the "0." version number.  I'm also still adding functionality,
so there may be interface changes from one version to the next if I decide I 
want to do something a different way.  Just be warned.  Any such changes will be 
mentioned in \S\ref{history} of subsequent releases.

The Template Matrix/Vector (TMV) Library is a C++ class library designed to make
writing code with vectors and matrices both transparent and fast.  Transparency 
means that when you look at your code months later, it is obvious what the code
does, making it easier to debug.  Fast means the execution time of the code should
be as fast as possible - this is mostly algorithm dependent, so we want the 
underlying library code to use the fastest algorithms possible.

If there were another free C++ Matrix library available that satisfied these requirements
and provided all (or even most) of the functionality I wanted, I probably would
not have written this.  But, at least when I started writing this, the available matrix libraries
were not very good.  Either they didn't have good operator overloading, or they 
didn't do complex matrices well, or they didn't include singular value decompositions,
or something.  Anyway, since I did decide to write 
my own library, hopefully other people can benefit from my efforts and will find 
this to be a useful product.

Given the above basic guidelines, the specific design features that I 
have incorporated into the code include:

\begin{enumerate}
\item
\textbf{Operator overloading}

Matrix equations look like real math equations in the code. 
For example, one can write 
\begin{verbatim}
v2 = m * v1;
v2 += m * v1;
m *= 3.;
v2 += m1*v1 + 3*v2 + m2.Transpose()*v3;
\end{verbatim} 
to perform the corresponding operations.

I also define division:
\begin{verbatim}
x = b/A;
\end{verbatim}
to mean solve the matrix equation $A x = b$.  If $A$ has more rows than columns,
then the solution will be a least-square solution.

\item
\textbf{Delayed evaluation}

Equations like the above 
\begin{verbatim}
v2 = m * v1;
v2 += m * v1;
\end{verbatim}
do not create a temporary vector before assigning or adding the result to \tt{v2}. 
However, a more complicated equation like
\begin{verbatim}
m2 += m1*v1 + 3*v2 + m2.Transpose()*v3;
\end{verbatim}
does not have a specialized routine, 
so it will require a couple temporary \tt{Vector}s.
Generally, if a statement performs just one operation, no temporary will be needed.  
Equations like this will give the right answer, but may not be quite as efficient as if you expand
the code to perform one operation per line.

\item
\textbf{Template classes}

Of course, all of our matrix and vector classes are templates, so we can have 
\begin{verbatim}
Matrix<float>
Matrix<double>
Matrix<complex<double> >
Matrix<long double>
Matrix<Quad>
\end{verbatim}
or whatever.

\item
\textbf{Mix complex/real}

One can multiply a real matrix by a complex
vector without having to copy the matrix to a new complex matrix.  Likewise
for the other arithmetic operations.  
However, we
do not allow mixing of underlying data types (float with double, for example), 
with the exception of simple assignments.

\item
\textbf{Views}

Operations like \tt{m.Tranpose()} or \tt{v.SubVector(3,8)}
return ``views'' of the underlying data rather than copying to new storage.  
This model helps delay calculations, which increases efficiency.  And the syntax
is fairly obvious.  For example:
\begin{verbatim}
v.SubVector(3,8) *= 3.;
m.row(3) += 4. * m.row(0);
m *= m.Transpose();
\end{verbatim}
modifies the underlying \tt{v} and \tt{m} in the obvious ways.

Note that in the last equation, \tt{m.Transpose()} uses the same storage as
\tt{m}, which is getting overwritten in the process.  The code recognizes 
this conflict and uses temporary storage to obtain the correct result.
However, if the overlapping storage is not at the first element, the code will not
recognize it, and you may get the wrong answer. 
See ``Alias checking'' below for more about this.

\item
\textbf{C- or Fortran-style indexing}

Both C- and Fortran-style (i.e. zero-based or one-based) indexing are possible for
element access of the matrices and vectors.

With C-style indexing, all matrix and vector indexing starts with 0.  
So the upper-left element of a matrix is \tt{m(0,0)}, not \tt{m(1,1)}. 
Likewise, the lower right element of an $M \times N$ matrix is \tt{m(M-1,N-1)}.
For element ranges, such as \tt{v.SubVector(0,10)}, the first number is the 
index of the first element, and the second number is ``one-past-the-end'' of 
the range.  So, this would return a 10 element vector from \tt{v(0)} to
\tt{v(9)} inclusive, not an 11 element one including \tt{v(11)}.

With Fortran-style indexing, all matrix and vector indexing starts with 1.
So the upper-left element of a matrix is \tt{m(1,1)}.  Likewise, the lower right
element of an $M \times N$ matrix is \tt{m(M,N)}.  For element ranges, such as
\tt{v.SubVector(1,10)}, the first number is the index of the first element, and
the second number is the last element.  So, this would return a 10 element vector
from \tt{v(1)} to \tt{v(10)} inclusive, which represents the same actual elements 
as the C-style example above.

\item
\textbf{Special matrices}

Many applications use matrices with known sparsity or a specific structure.  
The code is able to exploit a number of these structures for increased
efficiency in both speed and storage.  So far the following special matrices are
available: diagonal, upper/lower triangle, symmetric, hermitian, band, and 
symmetric or hermitian band.  
Special types of band matrices, such as 
upper and lower band, tridiagonal or Hessenberg, may all be declared as a 
regular band matrix.  The code checks the number of sub- and super-diagonals 
and uses the right algorithm when such a specialization is advantageous for a 
particular calculation.

\item
\textbf{Flexible storage}

Both row-major and column-major storage are possible as an optional
extra template parameter.
For band-diagonal matrices, we also allow diagonal-major storage.
This can aid I/O which may require a particular format to mesh with another
program.  Also, some algorithms 
are faster for one storage than than the other, so it can be worth switching storage
and checking the speed difference.

\item
\textbf{Alias checking}

Expressions such as \tt{m *= m} pose a problem for many matrix libraries, since
no matter what order you do the calculation, you will necessarily overwrite elements 
that you need for later stages of the calculation.  
The TMV code automatically recognizes the
conflict (generally known as an alias) and creates the needed temporary storage.

The code only checks the addresses of the first elements of the different objects.
So things like \tt{m = LowerTriMatrixViewOf(m) * UpperTriMatrixViewOf(m)} will work
correctly, even though there are three types of matrices involved, since the 
address of the upper-left corner of each matrix is the same.  (This particular
expression does not even need a temporary.  The code orders the steps of this
calculation so it can be done in place.)

However, \tt{v.SubVector(5,15) += v.SubVector(0,10)} will be calculated 
incorrectly, since the subvectors start at different locations, so the code
doesn't notice the aliasing.  Here, elements 5-9 will be overwritten before
they are added to the left-side vector.

Therefore, some care is still needed.  But this limited check is sufficient
for most applications.

\item
\textbf{BLAS}

For the combinations of types for which there are existing BLAS
routines, the code can call the optimized BLAS routines instead of its own 
code.  For other combinations (or for user defined types like 
\tt{Quad} or somesuch), 
the code does its best to order the steps of the calculation to be 
reasonably efficient, but it definitely is not as fast as BLAS for
most machines.

This feature can be turned off at compile time if desired with the 
compilation flag -DNOBLAS, although we do not generally recommend it if 
BLAS is available on your machine, and speed is important for your 
application.  For some operations, the BLAS routines can be a factor of 
10 or more faster than the native TMV code.  (I'm working on reducing this
gap, but I don't really expect to ever do better than a factor of 2 slower
than an optimized BLAS library.)

\item
\textbf{LAPACK}

When possible, the code can call LAPACK routines, which
may be faster than the native TMV code.  For types which don't have LAPACK routines, 
the code uses blocked and/or recursive algorithms which are similarly fast.  
Again, this feature can be turned off at compile time, this time with -DNOLAP
(-DNOBLAS necessarily turns off the LAPACK calls as well.)

For almost all algorithms, 
the TMV code is approximately as fast as LAPACK routines -
sometimes faster, since most
LAPACK distributions do not use recursive algorithms yet, which are generally
slightly faster on modern machines with good compilers. 
So if you don't want to deal with getting LAPACK
up and running, it won't generally be too bad, speedwise, 
to turn off the LAPACK calls.

Currently, the exception to this statement
is the singular value decomposition.  I have not yet 
implemented either the divide-and-conquer algorithm, nor the new RRR
(Relatively Robust Representation) algorithm by Dhillon.  So if the code is
spending a significant time doing SVD's, it may be worth having it call the
LAPACK routines.

\item
\textbf{Comments}

The code is highly commented, especially for the more complicated algorithms.
I have this egotistical quirk that I need to understand all of the code I write.
So none of the algorithms are just taken verbatim from LAPACK or anything like that.
I force myself to understand the algorithm (usually from Golub and van Loan,
but sometimes from a journal article or something I find on the net),
write a long comment explaining how it works, and then finally write the code
from scratch.  Then I usually compare my code to the LAPACK version, at which point
I have occasionally found ways to either speed it up or make it more accurate
(which also get commented, of course).

But the point is - if you want to understand how, say, the blocked Bunch-Kauffman
algorithm works, you'll probably do better looking at the comment in the 
TMV code than trying to decipher the LAPACK routines.  Plus, the code itself is
usually a lot more readable, since the arithmetic is written with the operator
overloads and such rather than the cryptic BLAS function names.

\end{enumerate}

All of the basic TMV classes and functions, including the \tt{Vector} and \tt{Matrix} 
classes, can be accessed with
\begin{verbatim}
#include "TMV.h"
\end{verbatim}
This file includes all the other files for the basic TMV routines.  Some of the the 
special matrices are not included in this.  See their sections below for the name 
of the file to include to access those classes.

All of the TMV classes and functions reside in the namespace \tt{tmv}. 
And of course, they are all templates.
So if you want to declare a 10x10 \tt{Matrix}, one would write:
\begin{verbatim}
tmv::Matrix<double> m(10,10);
\end{verbatim}

If writing \tt{tmv::} all the time is cumbersome, one can use \tt{using}
statements near the top of the code:
\begin{verbatim}
using tmv::Matrix;
using tmv::Vector;
\end{verbatim}
Or, while generally considered bad coding style, one can import the whole namespace:
\begin{verbatim}
using namespace tmv;
\end{verbatim}

Throughout most of the documentation, we will write \tt{T} for the underlying type.
Wherever you see \tt{T}, you should put \tt{double} or \tt{float} or whatever.
For a user defined type, like \tt{Quad}, for example, the main requirements are that 
both of the functions:
\begin{verbatim}
std::numeric_limits<T>::epsilon()
sqrt(T x) 
\end{verbatim}
be defined appropriately, where \tt{T} is your type name.  

Some functions will return a real value or require a real argument, even if \tt{T}
is complex.  In these cases, I will write \tt{RT} to indicate 
``the real type associated with \tt{T}''.

It may be worth noting that \tt{Matrix<int>} is possible as well.  
However, only simple 
routines like multiplication and addition will give correct answers.  If you try to 
divide by a \tt{Matrix<int>}, for example, the required calculations are impossible 
for \tt{int}'s,
so the result will not be correct.  But since the possibility of multiplication of 
integer Matrices seemed desirable, we do allow them to be used.  
\emph{Caveat programor}.  If debugging is turned on (or more accurately, not
turned off via the compile flag -DNDEBUG), then trying to do anything which requires
\tt{sqrt} or \tt{epsilon} for \tt{int}s will result in a runtime error.

\newpage
\section{Vectors}

The \tt{Vector} class is our mathematical vector.  Not to be confused with
the standard template library's \tt{vector} class.  
Our \tt{Vector} class name is capitalized, while the STL \tt{vector} is not.
If this is not enough of a difference for you, and you are using both extensively in your code,
we recommend keeping the full \tt{tmv::Vector} designation for ours to distinguish them.

\tt{Vector} inherits from its base class \tt{GenVector} (i.e. ``generic vector'').
Most operations which do not 
modify the data are actually defined in \tt{GenVector} rather than \tt{Vector}, although
some routines are overridden in \tt{Vector} for speed reasons.  

The other classes which inherit from \tt{GenVector} are \tt{VectorView}, 
\tt{ConstVectorView} (both described in more detail below - see \S\ref{VViews}), and
\tt{VectorComposite}, which is the base class for all arithmetic operations that
return a (logical) \tt{Vector}.  This means that any of these other objects can be used
any time a \tt{Vector} can be used in a non-assignable capacity.  For example,
\tt{Norm(v1+m*v2)} is completely valid, and will automatically create the necessary
temporary \tt{Vector} to store the result of the mathematical operation in the parentheses.

There is another template argument for a \tt{Vector} in addition to \tt{T} (which
represents the data type of the elements).  The second template argument
may be either \tt{tmv::CStyle} or \tt{tmv::FortranStyle}.  Or it may be omitted,
in which case \tt{CStyle} is assumed.  This argument governs how the element
access is performed.

With C-style indexing, the first element of a \tt{Vector} of length \tt{N} is 
\tt{v(0)} and the last is \tt{v(N-1)}.  Also, methods which take range arguments
use the common C convention of ``one-past-the-end'' for the last element;
so \tt{v.SubVector(0,3)} returns a 3-element vector, not 4.

With Fortran-style indexing, the first element of the vector is \tt{v(1)} and the 
last is \tt{v(N)}.  And ranges indicate the first and last elements, so the same
subvector as above would now be accessed using \tt{v.SubVector(1,3)} to return
the first three elements.

All views of a \tt{Vector} keep the same indexing style as the original unless you
explicitly change it with a cast.  You can cast a \tt{VectorView<T,CStyle>} as
a \tt{VectorView<T,FortranStyle>} and vice versa.  Likewise for
\tt{ConstVectorView}.  

The only thing to watch out for about the indexing is that \tt{GenVector}
and \tt{VectorComposite} do not have the extra \tt{index} 
template argument and are always
indexed using the C-style convention.  Therefore, 
if you want to index a \tt{GenVector}
using the Fortran-style convention, you would need to recast it
as an object of type
\tt{ConstVectorView<T,FortranStyle>}.  A \tt{VectorComposite} would 
not generally be indexed, but if you did want to do so using the 
Fortran-style conventions, you would need to explicitly instantiate it
as a \tt{Vector<T,FortranStyle>}.

\subsection{Constructors}

Here, \tt{T} is used to represent the data type of the elements of the \tt{Vector}
(e.g. \tt{double}, \tt{complex<double>}, \tt{int}, etc.) and \tt{index} is either
\tt{tmv::CStyle} or \tt{tmv::FortranStyle}.  In all of the constructors the
\tt{index} template argument may be omitted, in which case \tt{CStyle} is assumed.

\begin{itemize}
\item 
\begin{verbatim}
tmv::Vector<T,index> v(size_t n)
\end{verbatim}
Makes a \tt{Vector} of size \tt{n} with \underline{uninitialized} values.
If debugging is turned on (this is actually the default - 
turn off debugging by compiling
with -DNDEBUG), then the values are in fact initialized to 888.  
This should help you notice
when you have neglected to initialize the \tt{Vector} correctly.

\item
\begin{verbatim}
tmv::Vector<T,index> v(size_t n, T x)
\end{verbatim}
Makes a \tt{Vector} of size \tt{n} with all values equal to \tt{x}

\item
\begin{verbatim}
tmv::Vector<T,index> v(size_t n, const T* vv)
tmv::Vector<T,index> v(const std::vector<T>& vv)
\end{verbatim}
Makes a \tt{Vector} which copies the elements of \tt{vv}.
For the first one, \tt{n} specifies the length.  The other two get
the length from \tt{vv}.

\item
\begin{verbatim}
tmv::Vector<T,index> v = tmv::BasisVector<T,index>(size_t n, size_t i)
\end{verbatim}
Makes a \tt{Vector} whose elements are all \tt{0}, except \tt{v(i) = 1}.
Note the \tt{BasisVector} also has the \tt{index} template argument to 
indicate which element is meant by \tt{v(i)}.  Again, if it is omitted,
\tt{CStyle} is assumed.

\item
\begin{verbatim}
tmv::VectorView<T,index> v = 
        tmv::VectorViewOf(T* vv, size_t n, int step = 1)
tmv::ConstVectorView<T,index> v = 
        tmv::VectorViewOf(const T* vv, size_t n, int step = 1)
\end{verbatim}
Makes a \tt{VectorView} (see \S\ref{VViews} below) which refers to the exact
elements of \tt{vv}, not copying them to new storage.  The parameter \tt{n}
is the number of values to include in the view.
The optional \tt{step} parameter allows a non-unit spacing between 
successive vector elements in memory. 

\item
\begin{verbatim}
tmv::Vector<T,index> v1 = v2
tmv::Vector<T,index> v1(const tmv::GenVector<T2>& v2)
\end{verbatim}
Copy the \tt{Vector v2}, which may be of any type \tt{T2} so long
as values of type \tt{T2} are convertible into type \tt{T}.
The assignment operator has the same flexibility.

\end{itemize}

\subsection{Access}

\begin{itemize}
\item
\begin{verbatim}
v.size()
\end{verbatim}
Returns the size (length) of \tt{v}.

\item
\begin{verbatim}
v[i]
v(i)
\end{verbatim}
These are equivalent.  Each returns the \tt{i}-th element of \tt{v}.  
With \tt{index = CStyle}, the first 
element is \tt{v(0)}, and the last element is \tt{v(n-1)}.
With \tt{index = FortranStyle}, they are \tt{v(1)} and \tt{v(n)}.

If \tt{v} is a 
non-\tt{const Vector}, then the return type is a reference (\tt{T\&}).

If \tt{v} is a 
\tt{const Vector}, a \tt{ConstVectorView}, or a \tt{GenVector}, 
then the return type is just the value, not a reference.

If \tt{v} is a \tt{VectorView}, then the return type is an object which is
an lvalue (i.e. it is assignable), but which may not be \tt{T\&}.
Specifically, it has the type \tt{typename Vector<T>::reference}.
For a real-typed \tt{VectorView}, it is just \tt{T\&}.  But for a
complex-typed \tt{VectorView}, the return type is an object which keeps track of the
possibility of a conjugation.

\item
\begin{verbatim}
typename tmv::Vector<T>::iterator v.begin()
typename tmv::Vector<T>::iterator v.end()
typename tmv::Vector<T>::const_iterator v.begin() const
typename tmv::Vector<T>::const_iterator v.end() const
typename tmv::Vector<T>::reverse_iterator v.rbegin()
typename tmv::Vector<T>::reverse_iterator v.rend()
typename tmv::Vector<T>::const_reverse_iterator v.rbegin() const
typename tmv::Vector<T>::const_reverse_iterator v.rend() const
\end{verbatim}
These provide iterator-style access into a \tt{Vector}, which works just like
the standard template library's iterators.  If \tt{v} is a \tt{VectorView},
the iterator types are slightly different from the \tt{Vector} iterators, 
so you should declare them as:
\begin{verbatim}
typename tmv::VectorView<T>::iterator
\end{verbatim}
etc. instead.

\end{itemize}


\subsection{Views}
\label{VViews}

A \tt{VectorView<T>} object refers to the elements of some other object, such as a regular
\tt{Vector<T>} or \tt{Matrix<T>},
so that altering the elements in the view alters the
corresponding elements in the original object.  A \tt{VectorView}
can have non-unit
steps between elements (for example, a view of a column of a row-major
matrix).  It can also be a conjugation of the original
elements, so that
\begin{verbatim}
tmv::VectorView<double> cv = v.Conjugate();
cv(3) = z;
\end{verbatim}
would actually set the original element, \tt{v(3)} to \tt{conj(z)}. 

Also, we have to 
keep track of whether we are allowed to alter the original values or
just look at them.  
Since we want to be able to pass these views around, it turns out that
the usual \tt{const}-ing doesn't work the way you would want.
Thus, there are two objects which are views of a
\tt{Vector}:
\tt{ConstVectorView} and \tt{VectorView}.  
The first is only allowed to view,
not modify, the original elements.  The second is allowed to modify them.
This is akin to the \tt{const\_iterator} and \tt{iterator} types in the
standard template library.

One slightly non-intuitive thing about \tt{VectorView}s is that a 
\tt{const VectorView} is still mutable.  The \tt{const} in this case
means that one cannot change the components to which the view refers.
A \tt{VectorView} is inherently an object which can be used to 
modify the underlying data, regardless of any \tt{const} in front of it.

The following methods return views to portions of a \tt{Vector}.
If \tt{v} is either a (non-\tt{const}) \tt{Vector}
or a \tt{VectorView}, then a \tt{VectorView} is returned.
If \tt{v} is a \tt{const Vector}, a \tt{ConstVectorView}, or any other \tt{GenVector},
then a \tt{ConstVectorView} will be returned.  

\begin{itemize}
\item
\begin{verbatim}
v.SubVector(int i1, int i2, int istep=1)
\end{verbatim}
This returns a view to a subset of the original vector.
\tt{i1} is the first element in the subvector.
\tt{i2} is either ``one past the end'' (C-style) or the last element
(Fortran-style) of the subvector.
\tt{istep} is an optional step size.
Thus, if you have a \tt{Vector~v} of length 10, and you want to
multiply the first 3 elements by 2, with C-style indexing, you could write:
\begin{verbatim}
v.SubVector(0,3) *= 2.;
\end{verbatim}
To set all the even elements to 0, you could write:
\begin{verbatim}
v.SubVector(0,10,2).Zero();
\end{verbatim}
And then to output the last 4 elements of \tt{v}, you could write:
\begin{verbatim}
std::cout << v.SubVector(6,10);
\end{verbatim}

For Fortran-style indexing, the same steps would be accomplished by:
\begin{verbatim}
v.SubVector(1,3) *= 2.;
v.SubVector(1,9,2).Zero();
std::cout << v.SubVector(7,10);
\end{verbatim}

\item
\begin{verbatim}
v.Reverse()
\end{verbatim}
This returns a subvector view whose elements are the same as \tt{v},
but in the reverse order

\item
\begin{verbatim}
v.Conjugate()
\end{verbatim}
This returns the conjugate of a \tt{Vector} as a view, so it still points
to the same physical elements, but modifying this will set the 
actual elements in memory to the conjugate of what you set.  Likewise,
accessing an element will return the conjugate of the value in memory.

\item
\begin{verbatim}
v.View()
\end{verbatim}
Returns a view of a \tt{Vector}.  This seems like a silly function to have, but
if you write a function that takes a mutable \tt{Vector} argument, and you want to
be able to pass it views in addition to regular \tt{Vector}s, it is easier to write the 
function once with a \tt{VectorView} parameter.  Then you only need a second 
function with a \tt{Vector} parameter which calls the first function
using \tt{v.View()} as the argument:
\begin{verbatim}
double foo(const tmv::VectorView<double>& v)
{ ... [modifies v] ... }
double foo(tmv::Vector<double>& v)
{ return foo(v.View()); }
\end{verbatim}

If you are not going to be modifying \tt{v} in the function, you only need to
write one function, and you should use the base class \tt{GenVector} 
for the argument type:
\begin{verbatim}
double foo(const tmv::GenVector<double>& v)
{ ... [doesn't modify v] ... }
\end{verbatim}
The arguments could then be a \tt{const Vector}, a \tt{ConstVectorView},
or even a \tt{VectorComposite}.

\item
\begin{verbatim}
v.Real()
v.Imag()
\end{verbatim}
These return views to the real and imaginary parts of a complex \tt{Vector}.
Note the return type is a real view in each case:
\begin{verbatim}
tmv::Vector<std::complex<double> > v(10,std::complex<double>(1,4));
tmv::VectorView<double> vr = v.Real();
tmv::VectorView<double> vi = v.Imag();
\end{verbatim}

\end{itemize}

\subsection{Functions of a Vector }

Functions which do not modify the \tt{Vector} are defined in 
\tt{GenVector}, and so can be used for any type derived from \tt{GenVector}:
\tt{Vector, ConstVectorView, VectorView,} or \tt{VectorComposite}.

Functions which modify the \tt{Vector} are only defined for 
\tt{Vector} and \tt{VectorView}.

\subsubsection{Non-modifying functions}

Each of the following functions can be written in two ways, either as a method or a function.
For example, the expressions:
\begin{verbatim}
RT v.Norm()
RT Norm(v)
\end{verbatim}
are equivalent.  
In each case, \tt{v} can be any \tt{GenVector}.  Just to remind you, \tt{RT} refers to the 
real type associated with \tt{T}.  So for \tt{T} = \tt{double} or \tt{complex<double>},
\tt{RT} would be \tt{double}.

\begin{itemize}

\item
\begin{verbatim}
RT v.Norm1()
RT Norm1(v)
RT v.SumAbsElements()
RT SumAbsElements(v)
\end{verbatim}
The 1-norm of \tt{v}: $||v||_1 = \sum_i |v(i)|$.

\item
\begin{verbatim}
RT v.Norm2()
RT Norm2(v)
RT v.Norm()
RT Norm(v)
\end{verbatim}
The 2-norm of \tt{v}: $||v||_2 = (\sum_i |v(i)|^2)^{1/2}$.
This is the most common meaning for the norm of a vector, so we
define the \tt{Norm} function to be the same as \tt{Norm2}.

\item
\begin{verbatim}
RT v.NormSq()
RT NormSq(v)
\end{verbatim}
The square of the 2-norm of \tt{v}: $(||v||_2)^2 = \sum_i |v(i)|^2$.

\item
\begin{verbatim}
RT v.NormInf()
RT NormInf(v)
RT v.MaxAbsElement()
RT MaxAbsElement(v)
\end{verbatim}
The infinity-norm of \tt{v}: $||v||_\infty = \max_i |v(i)|$.

\item
\begin{verbatim}
RT v.MaxAbsElement(size_t* i=0);
RT MaxAbsElement(v, size_t* i=0);
RT v.MinAbsElement(size_t* i=0);
RT MinAbsElement(v, size_t* i=0);
T v.MaxElement(size_t* i=0);
T MaxElement(v, size_t* i=0);
T v.MinElement(size_t* i=0);
T MinElement(v, size_t* i=0);
\end{verbatim}
The maximum/minimum element either by absolute value (the first four cases) or actual value
(the last four cases).  For complex \tt{Vector}s, there is no way to define a 
max or min element, so just the real component of each element is used.
The \tt{i} argument is optional for all of these function.  
If it is present (and not 0), then \tt{*i} is set to the 
index of the max/min element returned.

\item
\begin{verbatim}
T v.SumElements()
T SumElements(v)
\end{verbatim}
The sum of the elements of \tt{v} $= \sum_i v(i)$.

\end{itemize}


\subsubsection{Modifying functions}

The following functions are methods of both \tt{Vector} and \tt{VectorView},
and they work the same way in the two cases, although there may be speed differences
between them.  

All of these are usually written on a line by themselves.  However, they do
return the (modified) \tt{Vector}, so you can string them together if you want.
For example:
\begin{verbatim}
v.Clip(1.e-10).ConjugateSelf().ReverseSelf();
\end{verbatim}
would first clip the elements at \tt{1.e-10}, then conjugate each element, then finally reverse the
order of the elements.  (This would probably not be considered 
very good programming style, however.)
Likewise, the expression:
\begin{verbatim}
foo(v.Clip(1.e-10));
\end{verbatim}
which would first clip the elements at \tt{1.e-10},
then pass the resulting \tt{Vector} to the function \tt{foo}.  

\begin{itemize}

\item
\begin{verbatim}
v.Zero();
\end{verbatim}
Clear the \tt{Vector v}.  i.e. Set each element to 0.

\item
\begin{verbatim}
v.SetAllTo(T x);
\end{verbatim}
Set each element to the value \tt{x}.

\item
\begin{verbatim}
v.Clip(RT thresh)
\end{verbatim}
Set each element whose absolute value is less than \tt{thresh} equal to 0.
Note that \tt{thresh} should be a real value even for complex valued
\tt{Vector}s.

\item
\begin{verbatim}
v.AddToAll(T x)
\end{verbatim}
Add the value \tt{x} to each element.

\item 
\begin{verbatim}
v.ConjugateSelf()
\end{verbatim}
Change each element into its complex conjugate.  
Note the difference between this and \tt{v.Conjugate()}, 
which returns a \underline{view} to a conjugated version of \tt{v} without
actually changing the underlying data.  This function, \tt{v.ConjugateSelf()},
does change the underlying data.

\item
\begin{verbatim}
v.ReverseSelf()
\end{verbatim}
Reverse the order of the elements.  Note the difference between this and 
\tt{v.Reverse()} which returns a \underline{view} to the elements in 
reversed order.

\item
\begin{verbatim}
v.MakeBasis(size_t i)
\end{verbatim}
Set all elements to 0, except for \tt{v(i)} = 1.

\item
\begin{verbatim}
v.Swap(size_t i1, size_t i2)
\end{verbatim}
Swap elements \tt{v(i1)} and \tt{v(i2)}.

\item
\begin{verbatim}
v.Permute(const size_t* p)
v.Permute(const size_t* p, size_t i1, size_t i2)
v.ReversePermute(const size_t* p)
v.ReversePermute(const size_t* p, size_t i1, size_t i2)
\end{verbatim}
The first one performs a series of swaps: (\tt{v(0)},\tt{v(p[0])}), (\tt{v(1)},\tt{v(p[1])}), ... 
The second starts at \tt{i1} and ends at \tt{i2-1} rather than
doing the whole range from 0 to \tt{n-1}.
The last two work the same way, but do the swaps in the opposite order.

Note: The indices listed in a permutation array (\tt{p}) always use 
the C-style convention, even if \tt{v} used Fortran-style indexing.

\item
\begin{verbatim}
v.Sort(size_t* p=0, tmv::ADType ad=tmv::ASCEND, 
        tmv::COMPType comp=tmv::REAL_COMP)
\end{verbatim}
Sorts the \tt{Vector~v}, returning the swaps required in the array \tt{p}.
If you do not care about the swaps, you can set \tt{p = 0}.
Note: the returned permutation array, \tt{p}, uses the C-style indexing convention
even if \tt{v} uses Fortran-style indexing.

The second parameter, \tt{ad}, determines whether the sorted \tt{Vector} 
will have its elements in ascending or descending order.  The possible values
are \tt{ASCEND} and \tt{DESCEND}. 

The third parameter, \tt{comp}, determines what component of the
elements to use for the sorting.  This is especially relevant if T is complex, 
since complex values are not intrinsically sortable.
The possible values are
\tt{REAL\_COMP}, \tt{ABS\_COMP},
\tt{IMAG\_COMP}, and \tt{ARG\_COMP}.
Only the first two make sense for non-complex \tt{Vector}s.  

\item
\begin{verbatim}
Swap(v1,v2)
\end{verbatim}
Swap the corresponding elements of \tt{v1} and \tt{v2}.  Note that this does physically
swap the data elements, not just some pointers to the data.  That's because this function
is mostly called on views into larger data fields: for example, swapping two rows of a 
\tt{Matrix}.  In any case, it always takes $O(N)$ time, never $O(1)$.

\end{itemize}

\subsection{Arithmetic}

\subsubsection{Operators}

All the usual operators work the way you would expect for \tt{Vector}s.  For shorthand in the 
following list, I use \tt{x} for a scalar of type \tt{T} or \tt{RT}, and
\tt{v} for a \tt{Vector}.  When there are two \tt{Vector}s
listed, they may either be both of the same type \tt{T}, or one may be of type 
\tt{T} and the other of \tt{complex<T>}.  Whenever \tt{v} is an lvalue,
if may be either a \tt{Vector} or a \tt{VectorView}.  Otherwise, it may be any \tt{GenVector}.

Also, I use the notation \tt{[+-]} to mean either \tt{+} or \tt{-}, since
the syntax is generally the same for these operators.
Likewise, I use \tt{[*/]} when their syntax is equivalent.

\begin{verbatim}
v2 = -v1;

v2 = x * v1;
v2 = v1 [*/] x;

v3 = v1 [+-] v2;

v [*/]= x;

v2 [+-]= v1;

x = v1 * v2;
\end{verbatim}
The last one, \tt{v1 * v2},
returns the inner product of two vectors, which is a scalar.  That is, the product
is a row vector times a column vector.  

This is the only case (so far)
where the specific row or column orientation of a vector matters.  For the others
listed here, the left side and the right side are implied to be of the same orientation, but
that orientation is otherwise arbitrary.
Later, when we get to a matrix times a vector, the orientation 
of the vector will be inferred from context.

\subsubsection{Subroutines}

Each of the above equations use deferred calculation so that the sum or product is not calculated
until the storage is known.  The equations can even be a bit more complicated without 
requiring a temporary.  Here are some equations that do not require a temporary 
\tt{Vector} for the calculation:

\begin{verbatim}
v2 = -(x1*v1 + x2*v2);
v2 += x1*(x2*(-v1));
v2 -= x1*(v1/=x2);
\end{verbatim}

The limit to how complicated the right hand side can be without using a 
temporary is set by the 
functions that the code eventually calls to perform the calculation.  
While you shouldn't ever
need to use these directly, it may help you understand when the code will require 
temporary \tt{Vector}s.  If you do use these, note that the \tt{v} parameters
are \tt{VectorView}s, rather than \tt{Vector}s.  So you would need to 
call them with \tt{v.View()} if \tt{v} is a \tt{Vector}.

\begin{itemize}

\item
\begin{verbatim}
MultXV(T x, const VectorView<T>& v)
\end{verbatim}
Performs the calculation \tt{v *= x}.

\item
\begin{verbatim}
MultXV(T x, const GenVector<T1>& v1, const VectorView<T>& v2)
\end{verbatim}
Performs the calculation \tt{v2 = x*v1}.

\item
\begin{verbatim}
AddVV(T x, const GenVector<T1>& v1, const VectorView<T>& v2)
\end{verbatim}
Performs the calculation \tt{v2 += x*v1}.

\item
\begin{verbatim}
AddVV(T x1, const GenVector<T1>& v1, T x2, const GenVector<T2>& v2,
        const VectorView<T>& v3)
\end{verbatim}
Performs the calculation \tt{v3 = x1*v1 + x2*v2}.

\item
\begin{verbatim}
T MultVV(const GenVector<T>& v1, const GenVector<T2>& v2)
\end{verbatim}
Performs the calculation \tt{v1*v2}.

\end{itemize}

More complicated arithmetic equations such as
\begin{verbatim}
v1 += x*(v1+v2+v3) + (x*v3-v1) 
\end{verbatim}
will require one or more temporary vectors, and so may be less efficient than 
you might like, but the code should return the correct result, no matter
how complicated the equation is.

\subsection{Input/Output}

The simplest output is the usual:
\begin{verbatim}
os << v
\end{verbatim}
where \tt{os} is any \tt{std::ostream}.
The output format is:
\begin{verbatim}
n ( v(0)  v(1)  v(2)  ...  v(3) )
\end{verbatim}
where \tt{n} is the length of the \tt{Vector}.

The same format can be read back in one of two ways:
\begin{verbatim}
tmv::Vector<T> v(n);
is >> v;
std::auto_ptr<tmv::Vector<T> > pv;
is >> pv;
\end{verbatim}
For the first version, the \tt{Vector} must already be declared, which 
requires knowing how big it needs to be.  If the input \tt{Vector} does not
match in size, an exception of type \tt{tmv::ReadError} is thrown.
The second version allows you to automatically get the size from the input.  
The \tt{Vector} pointed to by \tt{v2}
will be created with \tt{new} according to whatever size the input \tt{Vector} is.

Often, it is convenient to output only those values which aren't very small. 
This can be done using
\begin{verbatim}
v.Write(std::ostream& os, RT thresh)
\end{verbatim}
which is equivalent to
\begin{verbatim}
os << tmv::Vector<T>(v).Clip(thresh);
\end{verbatim}
but without requiring the temporary \tt{Vector}.

\newpage
\section{Dense Rectangular Matrices}

The \tt{Matrix} class is our dense matrix class.  It inherits from \tt{GenMatrix},
which (as for \tt{GenVector}) has the definitions of all the methods which 
do not modify the \tt{Matrix}.

The other classes which inherit from \tt{GenMatrix} are \tt{ConstMatrixView},
\tt{MatrixView} (see \S\ref{MViews} below), and \tt{MatrixComposite}, which 
is the base class for the various arithmetic operations which return a 
(logical) \tt{Matrix}.  

\tt{GenMatrix} in turn inherits from \tt{BaseMatrix}.  All of the various
special \tt{Matrix} classes also inherit from \tt{BaseMatrix}.
\tt{BaseMatrix} has virtual declarations for the functions which can be 
performed on any kind of \tt{Matrix} regardless of sparsity.

In addition to the data type template parameter (indicated here by \tt{T} as usual),
there is also a storage template parameter, which may be either 
\tt{RowMajor} or \tt{ColMajor}.

Finally, there is also a template argument indicating which indexing convention
you want the matrix to use, which may be either \tt{CStyle} or 
\tt{FortranStyle}.

With C-style indexing, the upper-left element of an $M \times N$ \tt{Matrix} is
\tt{m(0,0)}, the lower-left is \tt{m(M-1,0)}, the upper-right is \tt{m(0,N-1)},
and the lower-right is \tt{m(M-1,N-1)}.  Also, methods which take a pair of 
indices to define a range use the common C convention of ``one-past-the-end'' 
for the meaning of the second index.
So \tt{m.SubMatrix(0,3,3,6)} returns a $3 \times 3$ submatrix.

With Fortran-style indexing, the upper-left element of an $M \times N$ \tt{Matrix}
is \tt{m(1,1)}, the lower-left is \tt{m(M,1)}, the upper-right is \tt{m(1,N)},
and the lower-right is \tt{m(M,N)}.  Also, methods which take range arguments
take the pair of indices to be the actual first and last elements in the range.
So \tt{m.SubMatrix(1,3,4,6)} returns the same $3 \times 3$ submatrix as given above.

All views of a \tt{Matrix} keep the same indexing style as the original unless you
explicitly change it with a cast.  You can cast a \tt{MatrixView<T,CStyle>}
as a \tt{MatrixView<T,FortranStyle>} and vice versa.  (Likewise for 
\tt{ConstMatrixView}.)  However, as for \tt{GenVector}, you should be
aware that \tt{GenMatrix}
and \tt{MatrixComposite} do not have the extra template argument and are always
indexed using the C-style convention.  So if you want to index one of these 
using the Fortran-style convention, you need to (respectively) cast the 
\tt{GenMatrix} as a \tt{ConstMatrixView<T,FortranStyle>} or instantiate the
\tt{MatrixComposite} as a \tt{Matrix<T,FortranStyle>}.

You may omit the indexing template argument, in which case \tt{CStyle} is assumed.
And if so, you may also then omit the storage argument, in which case \tt{RowMajor}
is assumed.  If you want to specify \tt{FortranStyle} indexing, you need
to include the storage argument.

\subsection{Constructors}

We use \tt{stor} to indicate the storage template argument. 
This argument must be either 
\tt{tmv::RowMajor} or \tt{tmv::ColMajor}.  And \tt{index} indicates either 
\tt{tmv::CStyle} or \tt{tmv::FortranStyle}.  The default values for these
are \tt{tmv::RowMajor} and \tt{tmv::CStyle} if the arguments are omitted.

\begin{itemize}
\item 
\begin{verbatim}
tmv::Matrix<T,stor,index> m(size_t nrows, size_t ncols)
\end{verbatim}
Makes a \tt{Matrix} with \tt{nrows} rows and \tt{ncols} columns with 
\underline{uninitialized} values.
If debugging is turned on (this is actually the default -
turn off debugging by compiling
with -DNDEBUG), then the values are in fact initialized to 888.  
This should help you notice
when you have neglected to initialize the \tt{Matrix} correctly.

\item
\begin{verbatim}
tmv::Matrix<T,stor,index> m(size_t nrows, size_t ncols, T x)
\end{verbatim}
Makes a \tt{Matrix} with \tt{nrows} rows and \tt{ncols} columns with all 
values equal to \tt{x}

\item
\begin{verbatim}
tmv::Matrix<T,stor,index> m(size_t nrows, size_t ncols, const T* mm)
tmv::Matrix<T,stor,index> m(size_t nrows, size_t ncols, 
        const std::vector<T>& mm)
\end{verbatim}
Makes a \tt{Matrix} with \tt{nrows} rows and \tt{ncols} columns,
which copies the elements of \tt{mm}.

If \tt{stor} is \tt{RowMajor}, then the elements of \tt{mm} are taken to
be in row-major order: first the \tt{ncols} elements of the first row, then the
\tt{ncols} elements of the second row, and so on for the \tt{nrows} rows.
Likewise if \tt{stor} is \tt{ColMajor},
then the elements of \tt{mm} are taken to be in column-major order - 
the elements of the first column, then the second, and so on.

\item
\begin{verbatim}
tmv::Matrix<T,stor,index> m(const std::vector<std::vector<T> >& mm)
\end{verbatim}
Makes a \tt{Matrix} with elements \tt{m(i,j) = mm[i][j]}.  The size of
the \tt{Matrix} is taken from the sizes of the \tt{vector}s.  
(If \tt{index} is \tt{FortranStyle}, then \tt{m(i,j) = mm[i-1][j-1]}.)

\item
\begin{verbatim}
tmv::MatrixView<T,index> m = 
        tmv::MatrixViewOf(T* mm, size_t nrows, size_t ncols, 
        StorageType stor)
tmv::ConstMatrixView<T,index> m = 
        tmv::MatrixViewOf(const T* mm,  size_t nrows, size_t ncols, 
        StorageType stor)
\end{verbatim}
Makes a \tt{MatrixView} (see \S\ref{MViews} below) which refers to the exact
elements of \tt{mm}, not copying them to new storage.

\item
\begin{verbatim}
tmv::MatrixView<T,index> m = RowVectorViewOf(Vector<T>& v)
tmv::MatrixView<T,index> m = 
        RowVectorViewOf(const VectorView<T>& v)
tmv::ConstMatrixView<T,index> m = 
        RowVectorViewOf(const GenVector<T>& v)
\end{verbatim}
Makes a view of the \tt{Vector} which treats it as a $1\times n$ \tt{Matrix}
(i.e. a single row).  The first (non-const) version may also take a
\tt{VectorView} argument.

\item
\begin{verbatim}
tmv::MatrixView<T,index> m = ColVectorViewOf(Vector<T>& v)
tmv::MatrixView<T,index> m = 
        ColVectorViewOf(const VectorView<T>& v)
tmv::ConstMatrixView<T,index> m = 
        ColVectorViewOf(const Vector<T>& v)
\end{verbatim}
Makes a view of the \tt{Vector} which treats it as an $n \times 1$ \tt{Matrix}
(i.e. a single column).  The first (non-const) version may also take a
\tt{VectorView} argument.

\item
\begin{verbatim}
tmv::Matrix<T,index> m1 = m2
tmv::Matrix<T,index> m1(const GenMatrix<T2>& m2)
\end{verbatim}
Copy the \tt{Matrix m2}, which may be of any type \tt{T2} so long
as values of type \tt{T2} are convertible into type \tt{T}.
The assignment operator has the same flexibility.

\end{itemize}


\subsection{Access}

\begin{itemize}
\item
\begin{verbatim}
m.nrows() = m.colsize()
m.ncols() = m.rowsize()
\end{verbatim}
Returns the size of each dimension of \tt{m}. 

\item
\begin{verbatim}
m(i,j) = m[i][j]
\end{verbatim}
Returns the \tt{i,j} element of \tt{m}. i.e. the \tt{i}th element in the 
\tt{j}th column.  Or
equivalently, the \tt{j}th element in the \tt{i}th row. 

With C-style indexing, the upper-left element
of a \tt{Matrix} is \tt{m(0,0)}, the lower-left is \tt{m(nrows-1,0)},
The upper-right is \tt{m(0,ncols-1)}, and the lower-right is
\tt{m(nrows-1,ncols-1)}.

With Fortran-style indexing, these four elements would instead be
\tt{m(1,1)}, \tt{m(nrows,1)}, \tt{m(1,ncols)}, and \tt{m(nrows,ncols)},
respectively.

If \tt{m} is a 
non-\tt{const Matrix}, then the return type is a reference (\tt{T\&}).

If \tt{m} is a 
\tt{const Matrix}, a \tt{ConstMatrixView}, or a \tt{GenMatrix}, 
then the return type is just the value, not a reference.

If \tt{m} is a \tt{MatrixView}, then the return type is an object
which is an lvalue (i.e. it is assignable), but which may not be \tt{T\&}.
It has the type \tt{typename MatrixView<T>::reference}, 
(which is the same as \tt{typename VectorView<T>::reference}).
It is equal to \tt{T\&} for real \tt{MatrixView}s, but is
more complicated for complex \tt{MatrixView}s since it needs to 
keep track of the possibility of conjugation.

\item
\begin{verbatim}
m.row(size_t i)
m.col(size_t j)
\end{verbatim}
Return a view of the \tt{i}th row or \tt{j}th column respectively.
If \tt{m} is mutable (either a non-\tt{const Matrix} or a \tt{MatrixView}),
then a \tt{VectorView} is returned.  Otherwise, a \tt{ConstVectorView}
is returned.

\item
\begin{verbatim}
m.row(size_t i, size_t j1, size_t j2)
m.col(size_t j, size_t i1, size_t i2)
\end{verbatim}
Variations on the above, where only a portion of the row or column
is returned.  

For example, with C-style indexing, \tt{m.col(3,2,6)} returns a 4-element
vector view containing the elements
[\tt{m(2,3), m(3,3), m(4,3), m(5,3)}].

With Fortran-style indexing, the same elements are returned by \tt{m.col(4,3,6)}. 
(And these elements would be called: [\tt{m(3,4), m(4,4), m(5,4), m(6,4)}].)

\item
\begin{verbatim}
m.diag()
m.diag(int i)
m.diag(int i, size_t k1, size_t k2)
\end{verbatim}
Return the diagonal or one of the sub- or super-diagonals.
This first one returns the main diagonal.  For the second and third,
\tt{i=0} refers to the main diagonal; \tt{i>0} are the super-diagonals;
and \tt{i<0} are the sub-diagonals.  The last version returns a subset
of the diagonal with \tt{k1 = 0} or \tt{1} (for C-style or Fortran-style
respectively) starting at the left or top edge of 
the \tt{Matrix}.

\item
\begin{verbatim}
m.SubVector(size_t i, size_t j, int istep, int jstep, size_t size)
\end{verbatim}
If the above methods aren't sufficient to obtain the \tt{VectorView} you
need, we provide this
function which returns a view through the \tt{Matrix} 
starting at \tt{m(i,j)}, stepping by \tt{(istep,jstep)} between elements,
for a total of \tt{size} elements.  For example, the diagonal
from the lower-left to the upper-right of an $n \times n$ \tt{Matrix} 
would be obtained by: \tt{m.SubVector(n-1,0,-1,1,n)} for C-style or
\tt{m.SubVector(n,1,-1,1,n)} for Fortran-style.

\end{itemize}


\subsection{Views}
\label{MViews}

A \tt{MatrixView} object refers to some or all of the elements of a regular \tt{Matrix},
so that altering the elements in the view alters the
corresponding elements in the original object.  A \tt{MatrixView}
can be either row-major, column-major, or neither.  That is, the view can 
span a \tt{Matrix} with non-unit steps in both directions.
As with a \tt{VectorView}, it can also be a conjugation of the original
elements.

And again, just like a \tt{VectorView}, there are two view classes for
a \tt{Matrix}:
\tt{ConstMatrixView} and \tt{MatrixView}.  
The first is only allowed to view,
not modify, the original elements.  The second is allowed to modify them.

It is worth pointing out again that a 
\tt{const MatrixView} is still mutable, just like a \tt{const VectorView}.
The \tt{const} just means that you cannot change which elements
the view references.

The following methods return views to portions of a \tt{Matrix}.
If \tt{m} is either a (non-\tt{const}) \tt{Matrix}
or a \tt{MatrixView}, then a \tt{MatrixView} is returned.
If \tt{m} is a \tt{const Matrix}, \tt{ConstMatrixView}, or any other \tt{GenMatrix},
then a \tt{ConstMatrixView} will be returned. 

\begin{itemize}
\item
\begin{verbatim}
m.SubMatrix(int i1, int i2, int j1, int j2)
m.SubMatrix(int i1, int i2, int j1, int j2, int istep, int jstep)
\end{verbatim}
This returns a view to a submatrix contained within the original matrix.

If \tt{m} uses C-style indexing, 
the upper-left corner of the returned view is \tt{m(i1,j1)},
the lower-left corner is \tt{m(i2-1,j1)},
the upper-right corner is \tt{m(i1,j2-1)}, and
the lower-right corner is \tt{m(i2-1,j2-1)}.

If \tt{m} uses Fortran-style indexing, 
the upper-left corner of the view is \tt{m(i1,j1)},
the lower-left corner is \tt{m(i2,j1)},
the upper-right corner is \tt{m(i1,j2)}, and
the lower-right corner is \tt{m(i2,j2)}.

The second version allows for non-unit steps in the two directions.
To set a \tt{Matrix} to be a checkerboard of 1's, you could write 
(for C-style indexing):
\begin{verbatim}
tmv::Matrix<int> board(8,8,0)
board.SubMatrix(0,8,0,8,2,2).SetAllTo(1);
board.SubMatrix(1,9,1,9,2,2).SetAllTo(1);
\end{verbatim}

For Fortran-style indexing, the same thing would be accomplished by:
\begin{verbatim}
tmv::Matrix<int,tmv::RowMajor,tmv::FortranStyle> board(8,8,0)
board.SubMatrix(1,7,1,7,2,2).SetAllTo(1);
board.SubMatrix(2,8,2,8,2,2).SetAllTo(1);
\end{verbatim}

\item
\begin{verbatim}
m.Rows(size_t i1, size_t i2)
m.Cols(size_t j1, size_t j2)
\end{verbatim}
Since pulling out a bunch of contiguous rows or columns is a common 
submatrix use, we provide these functions.  They are shorthand for
\begin{verbatim}
m.SubMatrix(i1,i2,0,ncols)
m.SubMatrix(0,nrows,j1,j2)
\end{verbatim}
respectively.  (For Fortran-style indexing, replace the 0 with a 1.)

\item
\begin{verbatim}
m.RowPair(i1,i2)
m.ColPair(i1,i2)
\end{verbatim}
Another common submatrix is to select a pair of rows or columns, not 
necessarily adjacent to each other.  These are short hand for:
\begin{verbatim}
m.SubMatrix(i1,i2+(i2-i1),0,ncols,i2-i1,1)
m.SubMatrix(0,nrows,j1,j2+(j2-j1),1,j2-j1)
\end{verbatim}
respectively.  The equivalent in Fortran-style indexing would be:
\begin{verbatim}
m.SubMatrix(i1,i2,1,ncols,i2-i1,1)
m.SubMatrix(1,nrows,j1,j2,1,j2-j1).
\end{verbatim}

\item
\begin{verbatim}
m.Transpose()
m.Conjugate()
m.Adjoint()
\end{verbatim}
These return the transpose, conjugate, and adjoint (aka conjugate-transpose
or ``dagger'') of a \tt{Matrix}.  They point to the 
same physical elements as the original matrix, so modifying these will
correspondingly modify the original matrix.

Note that some people define the adjoint of a matrix as the determinant times
the inverse.  This combination is also called the adjugate or the cofactor matrix.
It is \underline{not} the same as our \tt{m.Adjoint()}.  Our \tt{Adjoint} is usually
written $m^\dagger$, and is sometimes referred to as the hermitian conjugate
or (rarely) tranjugate.  Our definition of the adjoint seems to be the more modern
usage.  Older texts tend to use the other definition.  However, if this is confusing
for you, it may be clearer to explicitly write out \tt{m.Conjugate().Transpose()},
which will not produce any efficiency reduction in your code over using
\tt{m.Adjoint()}.

\item
\begin{verbatim}
m.View()
\end{verbatim}
Returns a view of a \tt{Matrix}.  As with the \tt{Vector View()} function, it is mostly
useful for writing a function that can take either a \tt{Matrix\&} argument or a
\tt{const MatrixView\&} argument.  This lets you convert the first into the second.

\item
\begin{verbatim}
m.Real()
m.Imag()
\end{verbatim}
These return views to the real and imaginary parts of a complex \tt{Matrix}.
Note the return type is a real view in each case:
\begin{verbatim}
tmv::Matrix<std::complex<double> > m(10,std::complex<double>(1,4));
tmv::MatrixView<double> mr = m.Real();
tmv::MatrixView<double> mi = m.Imag();
\end{verbatim}

\end{itemize}

\subsection{Functions of a Matrix}

Functions which do not modify the \tt{Matrix} are defined in 
\tt{GenMatrix}, and so can be used for any type derived from \tt{GenMatrix}:
\tt{Matrix}, \tt{ConstMatrixView}, \tt{MatrixView}, or \tt{MatrixComposite}.
Functions which modify the \tt{Matrix} are only defined for 
\tt{Matrix} and \tt{MatrixView}.

\subsubsection{Non-modifying functions}

Each of the following functions can be written in two ways,
either as a method or a function.
For example, the expressions:
\begin{verbatim}
RT m.Norm()
RT Norm(m)
\end{verbatim}
are equivalent.  
In each case, \tt{m} can be any \tt{GenMatrix}.
As a reminder, \tt{RT} refers to the real type associated with \tt{T}.
In other words, \tt{T} is either the same as \tt{RT} or it is
\tt{std::complex<RT>}.

\begin{itemize}

\item
\begin{verbatim}
RT m.Norm1()
RT Norm1(m)
\end{verbatim}
The 1-norm of \tt{m}: 
$||m||_1 = \max_j (\sum_i |m(i,j)|)$.

\item
\begin{verbatim}
RT m.Norm2()
RT Norm2(m)
\end{verbatim}
The 2-norm of \tt{m}: $||m||_2 =$ the largest singular value of $m$, which
is also the square root of the largest eigenvalue of $(m^\dagger m)$
Be warned that this function can be fairly expensive if you have not
already performed an SV decomposition \tt{m}.

\item
\begin{verbatim}
RT m.NormInf()
RT NormInf(m)
\end{verbatim}
The infinity-norm of \tt{m}: 
$||m||_\infty = \max_i (\sum_j |m(i,j)|)$.

\item
\begin{verbatim}
RT m.NormF()
RT NormF(m)
RT m.Norm()
RT Norm(m)
\end{verbatim}
The Frobenius norm of \tt{m}: 
$||m||_F = (\sum_{i,j} |m(i,j)|^2)^{1/2}$.
This is the most common meaning for the norm of a matrix, so we
define the \tt{Norm} function to be the same as \tt{NormF}.

\item
\begin{verbatim}
RT m.NormSq()
RT NormSq(m)
\end{verbatim}
The square of the Frobenius norm of \tt{m}: 
$(||m||_F)^2 = \sum_{i,j} |m(i,j)|^2$.

\item
\begin{verbatim}
RT m.MaxAbsElement()
RT MaxAbsElement(m)
\end{verbatim}
The element of \tt{m} with the maximum absolute value: 
$||m||_\Delta = \max_{i,j} |m(i,j)|$.

\item
\begin{verbatim}
T m.Trace()
T Trace(m)
\end{verbatim}
The trace of \tt{m}: $\mathrm{Tr}(m) = \sum_i m(i,i)$.

\item
\begin{verbatim}
T m.Det()
T Det(m)
\end{verbatim}
The determinant of \tt{m}.  For speed issues regarding this function, see 
\S\ref{division} below on division.

\item
\begin{verbatim}
tmv::Matrix<T> minv = m.Inverse()
tmv::Matrix<T> minv = Inverse(m)
m.Inverse(Matrix<T>& minv)
\end{verbatim}
Set \tt{minv} to the inverse of \tt{m}.  
Again, see the division section for speed issues here.
If \tt{m} is not square, then \tt{minv} is set to the pseudo-inverse, or an approximate
pseudo-inverse.  If \tt{m} is singular, then an error may result, or the pseudo-inverse
may be returned, depending on the division method specified for the matrix.  
(See \S\ref{pseudoinverse} and \S\ref{singular} 
on pseudo-inverses and singular matrices below).

\item
\begin{verbatim}
m.InverseATA(Matrix<T>& cov);
\end{verbatim}
If \tt{m} has more rows than columns, then using it to solve a system of equations
really amounts to finding the least-square solution, since there is (typically) no
exact solution.  When you do this, \tt{m} is known as the ``design matrix'' of the system,
and is commonly called $A$.  Solving $A x=b$ gives $x$ as the least-square 
solution.  And the covariance matrix of the elements of $x$ is 
$\Sigma = (A^\dagger A)^{-1}$.
It turns out that computing this matrix is generally easy to do once you have 
performed the decomposition needed to solve for \tt{x} (either a QR or SV 
decomposition - see \S\ref{decompositions} below).  
Thus, we include this function which sets the
argument \tt{cov} to the equivalent of \tt{Inverse(m.Adjoint()*m)}, 
but generally does so much more efficiently than doing this explicitly, 
and also probably more accurately.

\end{itemize}

\subsubsection{Modifying functions}

The following functions are methods of both \tt{Matrix} and \tt{MatrixView},
and they work the same way for each.
As with the \tt{Vector} modifying functions, these all return a reference
to the newly modified \tt{Matrix}, so you can string them together if you want.

\begin{itemize}

\item
\begin{verbatim}
m.Zero();
\end{verbatim}
Clear the \tt{Matrix~m}.  i.e. Set each element to 0.

\item
\begin{verbatim}
m.SetAllTo(T x);
\end{verbatim}
Set each element to the value \tt{x}.

\item
\begin{verbatim}
m.Clip(RT thresh)
\end{verbatim}
Set each element whose absolute value is less than \tt{thresh} equal to 0.
Note that \tt{thresh} should be a real value even for complex valued
\tt{Matrix}es.

\item
\begin{verbatim}
m.SetToIdentity(T x = 1)
\end{verbatim}
Set to \tt{x} times the identity \tt{Matrix}.  
If the argument \tt{x} is omitted, it is 
taken to be \tt{1}, so \tt{m} is set to the identity \tt{Matrix}.
This is equivalent to \tt{m.Zero().diag().SetAllTo(x)}.

\item 
\begin{verbatim}
m.ConjugateSelf()
\end{verbatim}
Conjugate each element.  Note the difference between this and \tt{m.Conjugate()}, 
which returns a \underline{view} to the conjugate of \tt{m} without
actually changing the underlying data.  Contrariwise, \tt{m.ConjugateSelf()}
does change the underlying data.

\item
\begin{verbatim}
m.TransposeSelf()
\end{verbatim}
Transpose the \tt{Matrix}.  Note the difference between this and 
\tt{m.Transpose()} which returns a \underline{view} to the transpose without 
actually changing the underlying data.

\item
\begin{verbatim}
m.SwapRows(size_t i1, size_t i2)
m.SwapCols(size_t j1, size_t j2)
\end{verbatim}
Swap the corresponding elements of two rows or two columns.

\item
\begin{verbatim}
m.PermuteRows(const size_t* p)
m.PermuteRows(const size_t* p, size_t i1, size_t i2)
m.ReversePermuteRows(const size_t* p)
m.ReversePermuteRows(const size_t* p, size_t i1, size_t i2)
\end{verbatim}
These are equivalent to the corresponding routines for \tt{Vector}s 
(\tt{Permute} and \tt{ReversePermute}), performing a series of \tt{SwapRows} commands.

\item
\begin{verbatim}
m.PermuteCols(const size_t* p)
m.PermuteCols(const size_t* p, size_t j1, size_t j2)
m.ReversePermuteCols(const size_t* p)
m.ReversePermuteCols(const size_t* p, size_t j1, size_t j2)
\end{verbatim}
Same as above, but performing the a series of \tt{SwapCols} commands.

\item
\begin{verbatim}
Swap(m1,m2)
\end{verbatim}
Swap the corresponding elements of \tt{m1} and \tt{m2}.  Note that this does physically
swap the data elements, not just some pointers to the data, so it takes $O(N^2)$ time.

\end{itemize}

\subsection{Arithmetic}

\subsubsection{Operators}

We'll start with the simple operators that require little explanation aside from the 
notation: 
\tt{x} is a scalar, \tt{v} is a \tt{Vector}, and \tt{m} is a \tt{Matrix}:
\begin{verbatim}
m2 = -m1

m2 = x * m1
m2 = m1 [*/] x

m3 = m1 [+-] m2

m [*/]= x

m2 [+-]= m1

v2 = m * v1
v2 = v1 * m
v *= m

m3 = m1 * m2
m2 *= m1
\end{verbatim}
Note that the orientation of a vector is inferred from context.  \tt{m*v} involves a 
column vector, and \tt{v*m} involves a row vector.

For the product of two vectors, there are two orientation choices that make some sense.
It could mean a row vector times a column vector, which is called the inner product.
Or it could mean a column vector times a row vector, which is called the outer product.
We chose to let \tt{v1*v2} indicate the inner product, since this is far more common.
For the outer product, we use a different symbol:
\begin{verbatim}
m = v1 ^ v2
\end{verbatim}
One problem with this choice is that the order of operations for \tt{\^} 
is not the same as for \tt{*}.  
So, one needs to be careful when combining it with other calculations.
For example
\begin{verbatim}
m2 = m1 + v1 ^ v2   [ERROR!]
\end{verbatim}
will give a compile time error indicating that you can't add a 
\tt{Matrix} and a \tt{Vector}, 
because the operator \tt{+} has higher precedence for the compiler than \tt{\^}.
So you need to write:
\begin{verbatim}
m2 = m1 + (v1 ^ v2)
\end{verbatim}

Next, it is sometimes convenient to be able to treat scalars as the scalar times an
identity \tt{Matrix}.  So the following operations are allowed and use that convention:
\begin{verbatim}
m2 = m1 [+-] x
m2 = x [+-] m1
m [+-]= x
\end{verbatim}
For example, you could check if a \tt{Matrix} is numerically close to the identity matrix 
with:
\begin{verbatim}
if (Norm(m-1.) < 1.e-8) { [...] }
\end{verbatim}

\subsubsection{Subroutines}

As with \tt{Vector}s, we try to defer calculations until we have a place to store them.
So \tt{m*v} returns an object which can be assigned to a \tt{Vector},
but hasn't performed the calculation yet.  

The limit to how complicated the right hand side can be is set by the 
functions that the code eventually calls to perform the calculation.  
While you shouldn't ever
need to use these directly, it may help you understand when the code will create a
temporary \tt{Matrix}.  If you do use these, note that the last parameters are
\tt{MatrixView}s, rather than \tt{Matrix}es.  So you would need to 
call them with \tt{m.View()} if \tt{m} is a \tt{Matrix}.  (For \tt{MultMV},
it would be \tt{v.View()}.)

\begin{itemize}

\item
\begin{verbatim}
MultXM(T x, const MatrixView<T>& m)
\end{verbatim}
Performs the calculation \tt{m *= x}.

\item
\begin{verbatim}
MultMV<bool add>(T x, const GenMatrix<Tm>& m, 
        const GenVector<Tv1>& v1, const VectorView<T>& v2)
\end{verbatim}
Performs the calculation \tt{v2 (+=) x*m*v1} where ``(+=)'' means ``+='' 
if \tt{add} is true and ``='' if \tt{add} is false.

\item
\begin{verbatim}
AddMM(T x, const GenMatrix<T1>& m1, const MatrixView<T>& m2)
\end{verbatim}
Performs the calculation \tt{m2 += x*m1}.

\item
\begin{verbatim}
AddMM(T x1, const GenMatrix<T1>& m1, T x2, const GenMatrix<T2>& m2,
        const MatrixView<T>& m3)
\end{verbatim}
Performs the calculation \tt{m3 = x1*m1 + x2*m2}.

\item
\begin{verbatim}
MultMM<bool add>(T x, const GenMatrix<T1>& m1, 
        const GenMatrix<T2>& m2, const MatrixView<T>& m3)
\end{verbatim}
Performs the calculation \tt{m3 (+=) x*m1*m2} where ``(+=)'' means ``+='' 
if \tt{add} is true and ``='' if \tt{add} is false.

\item
\begin{verbatim}
Rank1Update<bool add>(T x, const GenVector<T1>& v1, 
        const GenVector<T2>& v2, const MatrixView<T>& m)
\end{verbatim}
Performs the calculation \tt{m (+=) x*(v1\^{}v2)} where ``(+=)'' means ``+='' 
if \tt{add} is true and ``='' if \tt{add} is false.

\end{itemize}

\subsection{Matrix Division}
\label{division}

One of the main things people want to do with a matrix is use it to solve a 
set of linear equations.  The set of equations can be written as a single matrix
equation:
\begin{equation}
\nonumber
A x = b
\end{equation}
where $A$ is a matrix and $x$ and $b$ are vectors.  $A$ and $b$ are known, 
and one want to solve for $x$.  Sometimes there are
multiple systems to be solved using the same coefficients, in which case
$x$ and $b$ become matrices as well.

\subsubsection{Operators}

Using the TMV classes, one would solve this equations by writing simply:
\begin{verbatim}
x = b / A
\end{verbatim}
Note that this really means $x = A^{-1} b$, which is different from $x = b A^{-1}$.
Writing the matrix equation as we did ($A x=b$) is much more common than 
swapping $A$ and $x$, so it makes sense to use this definition for the \tt{/}
operator.

However, we do allow for the possibility of wanting to right-multiply a vector
by $A^{-1}$ (in which case the vector is inferred to be a row-vector).  We designate
this operation by:
\begin{verbatim}
x = b % A
\end{verbatim}
which means $x = b A^{-1}$.

Given this explanation, the rest of the division operations should be self-explanatory,
where we use the notation \tt{[/\%]} to indicate that either \tt{/} or \tt{\%} may
be used with the above difference in meaning:
\begin{verbatim}
v2 = v1 [/%] m
m3 = m1 [/%] m2
v [/%]= m
m2 [/%]= m1
m2 = x [/%] m1
\end{verbatim}

If you feel uncomfortable using the \tt{/} and \tt{\%} symbols,
you can also explicitly write things like
\begin{verbatim}
v2 = m.Inverse() * v1
v3 = v1 * m.Inverse()
\end{verbatim}
which delay the calculation in exactly the same way that the above forms do.  
These forms
do not ever explicitly calculate the matrix inverse, since this is not (numerically) a
good way to perform these calculations.  Instead, the appropriate decomposition 
(see \S\ref{decompositions} below)
is used to calculate \tt{v2} and \tt{v3}.

\subsubsection{Least-square solutions}

If $A$ is not square, then the equation $A x = b$ does not have a unique solution.
If $A$ has more rows than columns (is ``tall''), then there is in general no solution.
And if $A$ has more columns than rows (is ``short''), then there are an infinite 
number of solutions.  

The tall case is more common.  In this case, one is not looking for an exact solution
for $x$, but rather the value of $x$ which minimizes $||b - A x||_2$.  This is the 
meaning of the least-square
solution, and is the value returned by \tt{x = b/A} for the TMV classes.

The short case is not so common, but can still be defined reasonably.  The
value returned by \tt{x = b/A} in this case is the value of \tt{x} which satisfies the 
equation and has minimum 2-norm, $||x||_2$.

When you have calculated a least-square solution for \tt{x}, it is common to want 
to know the covariance matrix of the returned values.  It turns out that this
matrix is $(A^\dagger A)^{-1}$.  It is not very efficient to calculate this matrix
explicitly and then invert it.  But once you have calculated the decomposition
needed for the division, it is quite easy.  So we provide the routine 
\begin{verbatim}
A.InverseATA(Matrix<T>& cov)
\end{verbatim}
to perform the calculation efficiently.  (Make sure you save the decomposition with 
\tt{A.SaveDiv()} - see \S\ref{efficiency} on efficiency below.)

\subsubsection{Decompositions}
\label{decompositions}

There are quite a few ways to go about solving the equations written above.  
The more efficient ways involve decomposing $A$ into a product of special 
matrices.  You can select which decomposition to use with the method:
\begin{verbatim}
m.DivideUsing(tmv::DivType dt)
\end{verbatim}
where \tt{dt} can be any of \{\tt{tmv::LU, tmv::QR, tmv::QRP, tmv::SV}\}.
If you do not specify which decomposition to use, \tt{LU} is the 
default for square matrices, and \tt{QR} is the default for
non-square matrices.

\begin{enumerate}
\item
\textbf{LU Decomposition}: 
(\tt{dt} = \tt{tmv::LU}) $A = P L U$, where $L$ is a lower-triangle 
matrix with all 1's along the diagonal, $U$ is an upper-triangle matrix, 
and $P$ is a permutation matrix.  This decomposition is only available for 
square matrices.

\item
\textbf{QR Decomposition}: 
(\tt{dt} = \tt{tmv::QR}) $A = Q R$ where $Q$ is a unitary matrix
and $R$ is an upper-triangle matrix.  (Note: real unitary matrices are
also known as orthogonal matrices.)  A unitary matrix is such that
$Q^\dagger Q = I$. 

If $A$ is tall with dimensions $M \times N$, 
then $R$ has dimensions $N \times N$, and $Q$
has dimensions $M \times N$.  In this case, $Q$ will
only be column-unitary.  That is $Q Q^\dagger \neq I$.

If $A$ is short, then $A^T$ is actually decomposed into $Q R$.

\item
\textbf{QRP Decomposition}: 
(\tt{dt} = \tt{tmv::QRP}) $A = Q R P$ where $Q$ is unitary, $R$ is 
upper-triangle, and $P$ is a permutation.  This decomposition is somewhat
slower than a simple QR decomposition, but it is numerically more stable if
$A$ is singular, or nearly so.  
(Singular matrices are discussed in more detail in \S\ref{singular} below.)

There are two slightly different algorithms for doing a QRP decomposition, controlled by a global
\tt{bool} variable: \tt{tmv::StrictQRP}.
If this is set to true, then the decomposition will make the diagonal elements
of $R$ be strictly decreasing (in absolute value) from upper-left to lower-right.

If it is false, however (the default), then there will be no diagonal element
of $R$ below and right of one which is \underline{much} smaller in absolute value,
where ``much'' means the ratio will be at most $\epsilon^{1/4}$, where 
$\epsilon$ is the machine precision for the type in question.

\item
\textbf{Singular Value Decomposition}: 
(\tt{dt} = \tt{tmv::SV}) $A = U S V$ where $U$ is unitary
(or column-unitary if $A$ is tall), $S$ is diagonal with all real, non-negative 
values, which decrease along the diagonal, and $V$ is unitary 
(or row-unitary if $A$ is short).  The singular value decomposition is most useful
for matrices which are singular or nearly so.  
We will discuss this decomposition in more detail in \S\ref{singular} on
singular matrices below.

\end{enumerate}

\subsubsection{Pseudo-inverse}
\label{pseudoinverse}

If \tt{m} is not square, then \tt{m.Inverse()} should return what is called the pseudo-inverse.
If \tt{m} is tall, then \tt{m.Inverse() * m} is the identity matrix.  But \tt{m * m.Inverse()}
is not an identity.  If \tt{m} is short, then the opposite holds.

Some features of the pseudo-inverse (here we use $X$ to represent the pseudo-inverse
of $M$):
\begin{align*}
M X M &= M \\
X M X &= X \\
(M X)^T &= M X \\
(X M)^T &= X M 
\end{align*}
For singular matrices (square or not),
one can also define a pseudo-inverse with the same properties.

In the first sentence of this section, I used the word ``should''.  
This is because the different decompositions calculate the pseudo-inverse 
differently and result in slightly
different answers.  For QR or QRP, the matrix returned by \tt{m.Inverse()} for
tall matrices isn't
quite correct.  It satisfies the first three of the above equations,
but not the last one
(for short M, the third equation fails).  For SV, however,
the returned matrix is the true pseudo-inverse and all four equations are satisfied.

For singular matrices, QR will fail to give a good pseudo-inverse, QRP will be close
(again failing only the last equation), and SV will be correct.

\subsubsection{Efficiency issues}
\label{efficiency}

Let's say you compute a matrix decomposition for a particular division calculation, and
then later want to use it again for a different right hand side:
\begin{verbatim}
x = b / m;
[...]
y = c / m;
\end{verbatim}
Ideally, the code would just use the same decomposition as had already been calculated
for \tt{x} when calculating \tt{y}, so this second division would be very fast.
However, what if somewhere in the \tt{[...]}, the \tt{Matrix m} was modified?  
Then using
the same decomposition would be incorrect - a new decomposition would be needed
for the new \tt{Matrix}.

One solution is to try to keep track of when a \tt{Matrix} gets changed.  
We could set an
internal flag whenever it is changed to indicate that
any existing decomposition is invalid.  While not impossible,
this type of check is made
quite difficult by the way different view objects can point to the same data.  
It would
be difficult to make sure that any change in one view invalidates the decompositions
of all other views to the same data.

Our solution is to err on the side of correctness over efficiency and
to always recalculate the decomposition by default.  Of course, this
can be quite inefficient, so we allow the programmer to override this
behavior for a specific \tt{Matrix} object with the method:
\begin{verbatim}
m.SaveDiv()
\end{verbatim}
After this call, whenever a decomposition is set, it is saved for any future uses.
You are telling the program to assume that the \tt{Matrix m} will not change anymore.

If you do modify \tt{m} after a call to \tt{SaveDiv()}, 
you can manually reset the decomposition with
\begin{verbatim}
m.ReSetDiv()
\end{verbatim}
which deletes any current saved decomposition, and recalculates it.
Similarly,
\begin{verbatim}
m.UnSetDiv()
\end{verbatim}
will delete any current saved decomposition, but not calculate a new one.
This can be used to free up the memory that the decomposition had been using.

Sometimes you may want to set a decomposition before you actually need it.
For example, the division may be in a speed critical part of the code, but you have
access to the \tt{Matrix} well before then.  You can tell the object to calculate 
the decomposition with
\begin{verbatim}
m.SetDiv()
\end{verbatim}
This may also be useful if you just want to perform and access the decomposition
separate from any actual division statement.

Also, if you change what kind of decomposition the \tt{Matrix} should use by calling
\tt{DivideUsing(...)}, then this will also delete any existing decomposition that 
might be saved.  (Unless you ``change'' it to the same thing, in which case it
leaves it there.)

Finally, there is another efficiency issue which may be important.  The default
behavior is to use extra memory for calculating the decomposition, so the original
\tt{Matrix} is left unchanged.  However, it is often the case that once you have 
calculated the decomposition, you don't need the original matrix anymore.
In that case, it is ok to overwrite the original matrix.  For very large matrices, 
the savings in memory may be significant.  (The $O(N^2)$ steps in copying 
the \tt{Matrix} is generally negligible compared to the $O(N^3)$ steps in performing
the decomposition.  So memory issues are probably the only reason to do this.)

Therefore, we provide another routine that lets the decomposition be 
calculated in place, overwriting the original \tt{Matrix}:
\begin{verbatim}
m.DivideInPlace()
\end{verbatim}

\subsubsection{Determinants}
\label{determinants}

The calculation of the determinant of a matrix requires similar calculations as
one of the above decompositions.  Since a determinant only really makes sense
for square matrices, one would typically perform an LU decomposition to calculate
the determinant.  Then the determinant of $A$ is just the determinant of $U$
(possibly times $-1$ depending on the details of $P$), which is simply the 
product of the values along the diagonal.

Therefore, calling \tt{m.Det()} involves calculating the LU decomposition, and then
finding the determinant of $U$.  
If you are also performing a division calculation, you should probably use
\tt{m.SaveDiv()} to avoid calculating the decomposition twice.

If you have set \tt{m} to use some other decomposition using 
\tt{m.DivideUsing(...)},
then the determinant will be determined from that decomposition instead 
(all of which are similarly easy).

\subsubsection{Accessing the decompositions}

Sometimes, you want to access the components of the decomposition directly,
rather than just use them for performing the division or calculating the determinant.

\begin{itemize}
\item
For the LU decomposition, we have:
\begin{verbatim}
ConstLowerTriMatrixView<T> m.LUD().GetL();
ConstUpperTriMatrixView<T> m.LUD().GetU();
size_t* m.LUD().GetP();
bool m.LUD().IsTrans();
\end{verbatim}
\tt{GetL()} and \tt{GetU()} return $L$ and $U$.
\tt{GetP()} returns the permutation array, $P$, in a form
which can be used as an argument 
to the routines like \tt{m.PermuteRows(P)}.  Finally, \tt{IsTrans()} returns whether
the decomposition is equal to \tt{m} or \tt{m.Transpose()}.  

The following should result in a \tt{Matrix m2} which is numerically very close to
the original \tt{Matrix m}:
\begin{verbatim}
tmv::Matrix<T> m2 = m.LUD().GetL() * m.LUD().GetU();
m2.ReversePermuteRows(m.LUD().GetP());
if (m.LUD().IsTrans()) m2.TransposeSelf();
\end{verbatim}

\item
For the QR decomposition, we have:
\begin{verbatim}
tmv::Matrix<T> m.QRD().GetQ();
tmv::ConstUpperTriMatrix<T>& m.QRD().GetR();
bool m.QRD().IsTrans();
\end{verbatim}
\tt{GetQ()} and \tt{GetR()} return $Q$ and $R$.  \tt{GetQ()} needs to make a new
\tt{Matrix}, since the \tt{QRD} class actually stores $Q$ in a compact form
rather than as the actual matrix.  This routine converts it into a regular matrix.  
(No such explicit conversion is required when the class is used to perform divisions.)
\tt{IsTrans()} returns whether
the decomposition is equal to \tt{m} or \tt{m.Transpose()}.  

The following should result in a \tt{Matrix m2} which is numerically very close to
the original \tt{Matrix m}:
\begin{verbatim}
tmv::Matrix<T> m2(m.nrows(),m.ncols());
tmv::MatrixView<T> m2v = 
        m.QRD().IsTrans() ? m2.Transpose() : m2.View();
m2v = m.QRD().GetQ() * m.QRD().GetR();
\end{verbatim}

\item
For the QRP decomposition, we have:
\begin{verbatim}
tmv::Matrix<T> m.QRPD().GetQ();
tmv::ConstUpperTriMatrixView<T>& m.QRPD().GetR();
size_t* m.QRPD().GetP();
bool m.QRPD().IsTrans();
\end{verbatim}
\tt{GetQ()} and \tt{GetR()} return $Q$ and $R$.
\tt{GetP()} returns the permutation array, $P$.
\tt{IsTrans()} returns whether
the decomposition is equal to \tt{m} or \tt{m.Transpose()}.  

The following should result in a \tt{Matrix m2} which is numerically very close to
the original \tt{Matrix m}:
\begin{verbatim}
tmv::Matrix<T> m2(m.nrows(),m.ncols());
tmv::MatrixView<T> m2v = 
        m.QRPD().IsTrans() ? m2.Transpose() : m2.View();
m2v = m.QRPD().GetQ() * m.QRPD().GetR();
m2v.ReversePermuteCols(m.QRPD().GetP())
\end{verbatim}

\item
For the SV decomposition, we have:
\begin{verbatim}
tmv::ConstMatrixView<T>& m.SVD().GetU();
tmv::ConstVectorView<RT>& m.SVD().GetS();
tmv::ConstMatrixView<T>& m.SVD().GetV();
bool m.SVD().IsTrans();
\end{verbatim}
\tt{GetU()} and \tt{GetV()} return $U$ and $V$.
\tt{GetS()} returns a vector, which is the diagonal
of the diagonal matrix $S$.  To use $S$ as a matrix, you can use
\tt{DiagMatrixViewOf(S)}.
\tt{IsTrans()} returns whether
the decomposition is equal to \tt{m} or \tt{m.Transpose()}.  

The following should result in a \tt{Matrix m2} which is numerically very close to
the original \tt{Matrix m}:
\begin{verbatim}
tmv::Matrix<T> m2(m.nrows(),m.ncols());
tmv::MatrixView<T> m2v = 
        m.SVD().IsTrans() ? m2.Transpose() : m2.View();
m2v = m.SVD().GetU() * DiagMatrixViewOf(m.SVD().GetS()) * 
        m.SVD().GetV();
\end{verbatim}
\end{itemize}

\subsubsection{Singular matrices}
\label{singular}

If a matrix is singular (i.e. its determinant is 0), 
then LU and QR decompositions
will fail when you attempt to divide by the matrix, since the 
calculation involves division by 0.  Numerically, a matrix
which is close to singular may end up with 0's from round-off
errors.  

Or, more commonly, a singular or nearly singular matrix
will just have numerically very small values rather than 0's.  In this 
case, there won't be an error, but the results will be numerically
very unstable, since the calculation will involve dividing by numbers
which are comparable to the machine precision, $\epsilon$.

You can check whether a Matrix is (exactly) singular with 
the method
\begin{verbatim}
bool m.Singular()
\end{verbatim}
which basically just returns whether \tt{m.Det() == 0}.  But this will not
tell you if a matrix is merely close to singular, so it does not guard 
against unreliable results.

Singular value decompositions provides a way to deal with 
nearly singular matrices.  There are a number of methods for the 
\tt{SVD} object which can help diagnose and fix 
potential problems with using a singular matrix.

First, the so-called ``singular values'', which are the elements of the
diagonal $S$ matrix, tell you how close the matrix is to being singular.
Specifically, if the ratio of the smallest and the largest singular values,
$S(N-1)/S(0)$, is close to the machine precision, $\epsilon$, for the 
underlying type (\tt{double}, \tt{float}, etc.),
then the matrix is singular, or nearly so.
Note: the value of $\epsilon$ is accessible with:
\begin{verbatim}
std::numeric_limits<T>::epsilon()
\end{verbatim}

The inverse of this ratio, $S(0)/S(N-1)$, is known as the ``condition'' of the matrix
(specifically the 2-condition, or $\kappa_2$), which can be obtained by:
\begin{verbatim}
m.SVD().Condition()
\end{verbatim}
The larger the condition, the closer the matrix is to singular, and the
less reliable any calculation would be.

So, how does the SVD help in this situation?  (So far we have diagnosed
the possible problem, but not fixed it.)

Well, we need to figure out what solution we want from a singular matrix.
If the matrix is singular, then there are an infinite number of solutions
to $A x = b$.  (Or, for a tall matrix, there are an infinite number of choices for 
$x$ which give the same minimum value of $||A x-b||_2$.)

Another way of looking at is is that there will be particular values of $y$ 
for which $A y = 0$.
Then given a solution $x$, the vector $x^\prime = x + \alpha y$ for any $\alpha$
will produce 
the same solution: $A x^\prime = A x$.

The usual desired solution is the $x$ with minimum 2-norm, $||x||_2$.
With the SVD, we can get this solution by setting to 0 all of the values
in $S^{-1}$ which would otherwise be infinity (or at least large compared
to $1/\epsilon$).  It is somewhat ironic that the best way to deal with 
an infinite value is to set it to 0, but that is actually the solution we want.

There are two methods which can be used to set which 
singular values to set to 0:
\begin{verbatim}
m.SVD().Thresh(RT thresh)
m.SVD().Top(size_t nsing)
\end{verbatim}
\tt{Thresh} sets to 0 any singular values with $S(i)/S(0) <$ \tt{thresh}.
\tt{Top} uses only the largest \tt{nsing} singular values, and sets the rest to 0.

The default behavior is equivalent to:
\begin{verbatim}
m.SVD().Thresh(std::numeric_limits<T>::epsilon());
\end{verbatim}
since at least these values are unreliable.  For different applications,
you may want to use a larger threshold value.

You can check how many singular values are currently considered non-zero with
\begin{verbatim}
size_t m.SVD().GetKMax()
\end{verbatim}

A QRP decomposition can deal with singular matrices similarly,
but it doesn't have the flexibility in checking for not-quite-singular,
but somewhat ill-conditioned matrices like the SVD does.
QRP will put all of the small elements of R's diagonal in the 
lower right corner.  Then it ignores any that are small compared to 
$\epsilon$ when doing the division.  For actually singular matrices,
this should produce essentially the same result as the SVD solution.

We should also mention again that the 2-norm of a matrix is the 
largest singular value, which is just $S(0)$ in our decomposition.
So this norm requires an $O(N^3)$ calculation to calculate the 
SVD, rather than the $O(N^2)$ calculation that the other norms need.
This is fairly expensive if you have not already calculated the SVD
(and saved it with \tt{m.SaveDiv()}).

If you just want to calculate the 2-norm, or even all of the singular values,
but don't need to do the actual division, then you don't need to 
accumulate the $U$ and $V$ matrices.  This saves a lot of the 
calculation time.  Or you may want $U$ or $V$, but not both for 
some purpose.  We provide some other \tt{DivType}s for these cases:
\begin{verbatim}
m.DivideUsing(tmv::SVS)
m.DivideUsing(tmv::SVU)
m.DivideUsing(tmv::SVV)
\end{verbatim}
\tt{SVS} only keeps $S$.  
\tt{SVU} keeps $S$ and $U$, but not $V$.
\tt{SVV} keeps $S$ and $V$, but not $U$.
For all three, an error will result if you try to then use the decomposition 
for division.

\subsection{Input/Output}

The simplest output is the usual:
\begin{verbatim}
os << m
\end{verbatim}
where \tt{os} is any \tt{std::ostream}.
The output format is:
\begin{verbatim}
nrows ncols
( m(0,0)  m(0,1)  m(0,2)  ...  m(0,ncols-1) )
( m(1,0)  m(1,1)  m(1,2)  ...  m(1,ncols-1) )
...
( m(nrows-1,0) ...  ...  m(nrows-1,ncols-1) )
\end{verbatim}

The same format can be read back in one of two ways:
\begin{verbatim}
tmv::Matrix<T> m(nrows,ncols);
is >> m;
std::auto_ptr<tmv::Matrix<T> > pm;
is >> pm;
\end{verbatim}
For the first version, the \tt{Matrix} must already be declared, which 
requires knowing how big it needs to be.  If the input \tt{Matrix} does not
match in size, a runtime error will occur.
The second version allows you to get the size from the input.  \tt{m2}
will be created (with \tt{new}) 
according to whatever size the input \tt{Matrix} is.

Often, it is convenient to output only those values which aren't very small. 
This can be done using
\begin{verbatim}
m.Write(std::ostream& os, RT thresh)
\end{verbatim}
which is equivalent to
\begin{verbatim}
os << tmv::Matrix<T>(m).Clip(thresh);
\end{verbatim}
but without requiring the temporary \tt{Matrix}.

\subsection{Small Matrices}

The algorithms for regular \tt{Matrix} operations are optimized to be
fast for large matrices.  Usually, this makes sense, since any code with
both large and small matrices will probably have its performance dominated
by the speed of the large matrix algorithms.

However, it may be the case that a particular program spends all of its
time using $2 \times 2$ or $3 \times 3$ matrices.  In this case, 
many of the features of the TMV code are undesirable.  For example, 
the alias checking in the operation \tt{v2 = m * v1} becomes a significant
fraction of the operating time.  Even the simple act of performing a function
call, rather than doing the calculation inline may be a big performance hit.

So we include the alternate matrix class called \tt{SmallMatrix}, along 
with the corresponding vector class called \tt{SmallVector}, for which most of 
the operations are done inline.  Furthermore, the sizes are template arguments,
rather than normal parameters.  This allows the compiler to easily optimize
simple calculations that might only be 2 or 3 arithmetic operations, which
may significantly speed up your code.

The class \tt{SmallMatrix} inherits from \tt{GenSmallMatrix},
which in turn inherits from \tt{BaseMatrix}.
Likewise, the class \tt{SmallVector} inherits from \tt{GenSmallVector}.

All the \tt{SmallMatrix} and \tt{SmallVector} routines are included by:
\begin{verbatim}
#include "TMV_Small.h"
\end{verbatim}

\subsubsection{Constructors}

\begin{itemize}
\item 
\begin{verbatim}
tmv::SmallVector<T,N,index> v()
\end{verbatim}
Makes a \tt{SmallVector} of size \tt{N}
with \underline{uninitialized} values.
If debugging is turned on (i.e. not turned off
with -DNDEBUG), then the values are in fact initialized to 888. 

\item
\begin{verbatim}
tmv::SmallVector<T,N,index> v(T x)
\end{verbatim}
Makes a \tt{SmallVector} with all values equal to \tt{x}.

\item
\begin{verbatim}
tmv::SmallVector<T,N,index> v(const T* vv)
tmv::SmallVector<T,N,index> v(const std::vector<T>& vv)
\end{verbatim}
Makes a \tt{SmallVector} with values copied from vv.

\item 
\begin{verbatim}
tmv::SmallVector<T,M,N,stor,index> v(const GenVector<T>& vv)
\end{verbatim}
Makes a \tt{SmallVector} from a regular \tt{Vector}.

\item
\begin{verbatim}
tmv::SmallVector<T,M,N,stor,index> v1 = v2
tmv::SmallVector<T,M,N,stor,index> v1(const GenSmallVector<T2>& v2)
\end{verbatim}
Copy the \tt{SmallVector v2}, which may be of any type \tt{T2} so long
as values of type \tt{T2} are convertible into type \tt{T}.
The assignment operator has the same flexibility.

\item 
\begin{verbatim}
tmv::SmallMatrix<T,M,N,stor,index> m()
\end{verbatim}
Makes an \tt{M} $\times$ \tt{N} \tt{SmallMatrix} 
with \underline{uninitialized} values.
If debugging is turned on (i.e. not turned off
with -DNDEBUG), then the values are in fact initialized to 888. 

\item
\begin{verbatim}
tmv::SmallMatrix<T,M,N,stor,index> m(T x)
\end{verbatim}
Makes an \tt{M} $\times$ \tt{N} \tt{SmallMatrix} with all values equal to \tt{x}.

\item
\begin{verbatim}
tmv::SmallMatrix<T,M,N,stor,index> m(const T* mm)
tmv::SmallMatrix<T,M,N,stor,index> m(const std::vector<T>& mm)
tmv::SmallMatrix<T,M,N,stor,index> m(
        const std::vector<std::vector<T> >& mm)
\end{verbatim}
Makes an \tt{M} $\times$ \tt{N} \tt{SmallMatrix} with values copied from mm.

\item 
\begin{verbatim}
tmv::SmallMatrix<T,M,N,stor,index> m(const GenMatrix<T>& mm)
\end{verbatim}
Makes a \tt{SmallMatrix} from a regular \tt{Matrix}.

\item
\begin{verbatim}
tmv::SmallMatrix<T,M,N,stor,index> m1 = m2
tmv::SmallMatrix<T,M,N,stor,index> m1(const GenSmallMatrix<T2>& m2)
\end{verbatim}
Copy the \tt{SmallMatrix m2}, which may be of any type \tt{T2} so long
as values of type \tt{T2} are convertible into type \tt{T}.
The assignment operator has the same flexibility.

\item
\begin{verbatim}
tmv::ConstSmallVectorView<T,N,step,isconj,index> v(
        const ConstVectorView<T>& v2)
tmv::SmallVectorView<T,N,step,isconj,index> v(
        const VectorView<T>& v2)
tmv::ConstSmallVectorView<T,M,N,stepi,stepj,isconj,index> m(
        const ConstMatrixView<T>& m2)
tmv::SmallVectorView<T,M,N,stepi,stepj,isconj,index> m(
        const MatrixView<T>& m2)
\end{verbatim}
Convert a regular \tt{Vector} or \tt{Matrix} view into a
small variety.  Note that the step size between elements along a
column (\tt{stepi}) and along a row (\tt{stepj}) must be
specified in the template arguments.  Also whether or not the
view is the conjugate of the actual data values (\tt{isconj}).
The parameters \tt{isconj} and \tt{index} may be omitted,
with the defaults \tt{false} and \tt{CStyle} being assumed.

\item
\begin{verbatim}
tmv::ConstSmallVectorView<T,N,1,false,I> v =
        tmv::SmallVectorViewOf<N,I>(const T* v);
tmv::SmallVectorView<T,N,1,false,I> v =
        tmv::SmallVectorViewOf<N,I>(T* v);
tmv::ConstSmallMatrixView<T,M,N,N,1,false,I> v =
        tmv::SmallVectorViewOf<M,N,RowMajor,I>(const T* v);
tmv::ConstSmallMatrixView<T,M,N,1,M,false,I> v =
        tmv::SmallVectorViewOf<M,N,ColMajor,I>(const T* v);
tmv::SmallMatrixView<T,M,N,N,1,false,I> v =
        tmv::SmallVectorViewOf<M,N,RowMajor,I>(T* v);
tmv::SmallMatrixView<T,M,N,1,M,false,I> v =
        tmv::SmallVectorViewOf<M,N,ColMajor,I>(T* v);
\end{verbatim}
View actual data as a \tt{SmallVector} or \tt{SmallMatrix}.
In all cases, the \tt{I} template parameter may be omitted, with
\tt{CStyle} being assumed.

\end{itemize}


\subsubsection{Access}

\begin{verbatim}
v.size()
v[i] = v(i)
m.nrows() = m.ncols() = m.colsize() = m.rowsize()
m(i,j) = m[i][j]
m.row(i)
m.col(j)
m.diag()
\end{verbatim}
The last three all return \tt{SmallVectorView}s.

\begin{verbatim}
v.SubVector(i1,i2)
m.row(i,j1,j2)
m.col(j,i1,i2)
m.diag(i)
m.diag(i,j1,j2)
m.SubVector(i,j,istep,jstep,size)
m.SubMatrix(i1,i2,j1,j2)
m.SubMatrix(i1,i2,j1,j2,istep,jstep)
m.Rows(i1,i2)
m.Cols(j1,j2)
m.RowPair(i1,i2)
m.ColPair(j1,j2)
\end{verbatim}
Since the sizes of all of these views are not known at compile time, 
these methods return a regular \tt{VectorView} or \tt{MatrixView} as
appropriate.  If you want, you can convert these into 
either a \tt{SmallVectorView}
or a \tt{SmallMatrixView} with an explicit cast as described above.

\begin{verbatim}
m.Transpose()
m.Conjugate()
m.Adjoint()
m.View()
m.Real()
m.Imag()
\end{verbatim}
These all return \tt{SmallMatrixView}s.

\subsubsection{Functions}

\tt{SmallVector}s and \tt{SmallMatrix}es all have exactly the same
function methods as the regular varieties.  Likewise, the syntax of the
arithmetic is identical.  There are only a few methods which 
are not done inline.  

First, I have not implemented any division routines inline yet.  So 
all divisions by a \tt{SmallMatrix} just use the regular 
\tt{Matrix} division routines.  This includes the related 
routines like \tt{Det} and \tt{Norm2}.

Second, reading a \tt{SmallMatrix} or \tt{SmallVector} from a file 
uses the regular
\tt{Matrix} I/O methods.  Also, there is no \tt{auto\_ptr} 
version of these read operations, since you need to know the 
size of a \tt{SmallMatrix} or \tt{SmallVector} at compile time anyway, 
so there is no way to wait until the file is read to determine the size.

Last, the \tt{Sort} command for a \tt{SmallVector} just uses the
regular \tt{Vector} version, since it calls the Standard Template
Library's \tt{sort} function anyway, so I don't think there would 
be any real advantage to inlining it.

\newpage
\section{Special Matrices}

Most functions and methods for the various special matrix varieties work the same
as the for regular dense, rectangular matrices.  
In these cases, we will just list the functions
that are allowed for each special matrix, with the implicit assumption that the 
effect is the same as for a regular \tt{Matrix}.  Of course, there are usually 
algorithmic speed-ups which the code will use to take advantage of the 
particular structure.
Whenever there is a difference in how a function works,
we will explain the difference.

\subsection{Diagonal matrices}

The \tt{DiagMatrix} class is our diagonal matrix class.  
A diagonal matrix is only non-zero
along the main diagonal of the matrix.  

The class \tt{DiagMatrix} inherits from \tt{GenDiagMatrix},
which in turn inherits from \tt{BaseMatrix}.
The various views and composite classes described below
also inherit from \tt{GenDiagMatrix}.

All the \tt{DiagMatrix} routines are included by:
\begin{verbatim}
#include "TMV_Diag.h"
\end{verbatim}

\subsubsection{Constructors}

As usual, the optional \tt{index} template argument specifies which indexing
style to use.

\begin{itemize}
\item 
\begin{verbatim}
tmv::DiagMatrix<T,index> m(size_t n)
\end{verbatim}
Makes an \tt{n} $\times$ \tt{n} \tt{DiagMatrix} with \underline{uninitialized} values.
If debugging is turned on (i.e. not turned off
with -DNDEBUG), then the values are in fact initialized to 888. 

\item
\begin{verbatim}
tmv::DiagMatrix<T,index> m(size_t n, T x)
\end{verbatim}
Makes an \tt{n} $\times$ \tt{n} \tt{DiagMatrix} with all values equal to \tt{x}.

\item
\begin{verbatim}
tmv::DiagMatrix<T,index> m(const GenVector<T>& v)
\end{verbatim}
Makes a \tt{DiagMatrix} with \tt{v} as the diagonal.

\item 
\begin{verbatim}
tmv::DiagMatrix<T,index> m(const GenMatrix<T>& mm)
\end{verbatim}
Makes a \tt{DiagMatrix} with the diagonal of \tt{mm} as the diagonal.

\item
\begin{verbatim}
tmv::DiagMatrixView<T,index> m = DiagMatrixViewOf(Matrix<T>& mm)
tmv::DiagMatrixView<T,index> m = 
        DiagMatrixViewOf(const MatrixView<T>& mm)
tmv::ConstDiagMatrixView<T,index> m = 
        DiagMatrixViewOf(const GenMatrix<T>& mm)
\end{verbatim}
Makes a \tt{DiagMatrixView} of the diagonal portion of \tt{mm} (namely the main 
diagonal).

\item
\begin{verbatim}
tmv::DiagMatrixView<T,index> m = DiagMatrixViewOf(Vector<T>& v)
tmv::DiagMatrixView<T,index> m = 
        DiagMatrixViewOf(const VectorView<T>& v)
tmv::ConstDiagMatrixView<T,index> m = 
        DiagMatrixViewOf(const GenVector<T>& v)
\end{verbatim}
Makes a \tt{DiagMatrixView} with \tt{v} as the diagonal of the matrix.

\item
\begin{verbatim}
tmv::DiagMatrixView<T,index> m = 
        tmv::DiagMatrixViewOf(T* vv, size_t n)
tmv::ConstDiagMatrixView<T,index> m = 
        tmv::DiagMatrixViewOf(const T* vv, size_t n)
\end{verbatim}
Make a \tt{DiagMatrix} with the actual memory elements, \tt{vv} as the diagonal.

\item
\begin{verbatim}
tmv::DiagMatrix<T,index> m1 = m2
tmv::DiagMatrix<T,index> m1(const GenDiagMatrix<T2>& m2)
\end{verbatim}
Copy the \tt{DiagMatrix m2}, which may be of any type \tt{T2} so long
as values of type \tt{T2} are convertible into type \tt{T}.
The assignment operator has the same flexibility.

\end{itemize}


\subsubsection{Access}

\begin{verbatim}
m.nrows() = m.ncols() = m.colsize() = m.rowsize() = m.size()
m(i,j)
m(i) = m(i,i)
\end{verbatim}
For the mutable \tt{m(i,j)} version, 
\tt{i} must equal \tt{j}.
If \tt{m} is not mutable, then \tt{m(i,j)} with \tt{i}$\neq$\tt{j} returns the 
value 0.

\begin{verbatim}
m.diag()
\end{verbatim}

\begin{verbatim}
m.SubDiagMatrix(int i1, int i2, int istep = 1)
\end{verbatim}
This is equivalent to \tt{DiagMatrixViewOf(m.diag().SubVector(i,j,istep))}.
\begin{verbatim}
m.Real()
m.Imag()
m.Transpose()
m.Conjugate()
m.Adjoint()
m.View()
\end{verbatim}

\subsubsection{Functions}

\begin{verbatim}
RT m.Norm1() = Norm1(m)
RT m.Norm2() = Norm2(m)
RT m.NormInf() = NormInf(m)
RT m.NormF() = NormF(m) = m.Norm() = Norm(m)
RT m.MaxAbsElement() = MaxAbsElement(m)
\end{verbatim}
(Actually for a diagonal matrix, all of these norms are equal.)
\begin{verbatim}
RT m.NormSq() = NormSq(m)
T m.Trace() = Trace(m)
T m.Det() = Det(m)
minv = m.Inverse() = Inverse(m)
m.Inverse(Matrix<T>& minv)
m.Inverse(DiagMatrix<T>& minv)
m.InverseATA(Matrix<T>& cov)
m.InverseATA(DiagMatrix<T>& cov)
\end{verbatim}
Since the inverse of a \tt{DiagMatrix} is a \tt{DiagMatrix},
we also provide a version of the \tt{Inverse} syntax which allows minv
to be a \tt{DiagMatrix}.  Likewise for \tt{InverseATA}.

\begin{verbatim}
m.Zero()
m.SetAllTo(T x)
m.Clip(RT thresh)
m.SetToIdentity(T x = 1)
m.ConjugateSelf()
m.TransposeSelf()
Swap(m1,m2)
\end{verbatim}

\subsubsection{Arithmetic}

In addition to \tt{x}, \tt{v}, and \tt{m} from before, we now add \tt{d} for a \tt{DiagMatrix}.

\begin{verbatim}
d2 = -d1

d2 = x * d1
d2 = d1 [*/] x

d3 = d1 [+-] d2
m2 = m1 [+-] d
m2 = d [+-] m1

d [*/]= x

d2 [+-]= d1
m [+-]= d

v2 = d * v1
v2 = v1 * d
v *= d

d3 = d1 * d2
m2 = d * m1
m2 = m1 * d
d2 *= d1
m *= d

d2 = d1 [+-] x
d2 = x [+-] d1
d [+-]= x
\end{verbatim}

\subsubsection{Division}

The division operations are:
\begin{verbatim}
v2 = v1 [/%] d
m2 = m1 [/%] d
m2 = d [/%] m1
d3 = d1 [/%] d2
v [/%]= d
d2 [/%]= d1
m [/%]= d
\end{verbatim}

There is only one allowed \tt{DivType} for a \tt{DiagMatrix}: \tt{LU}. 
And, since it is also the default behavior,
there is no reason to ever use this function.  Furthermore, since a \tt{DiagMatrix}
is already a $U$ matrix, the decomposition requires no work at all.
Hence, it is always done in place; no extra storage is needed, and 
the methods \tt{m.DivideInPlace()}, \tt{m.SaveDiv()}, etc. are irrelevant.

If a \tt{DiagMatrix} is singular, you can check easily with \tt{m.Singular()},
but there is no direct way to treat the division as a SVD and skip any
divisions by 0.  I suppose it would be easy to add that functionality, since a 
\tt{DiagMatrix} is also already an SV decomposition, ignoring the fact that
$S$ should be all real and positive, which is a small detail.  But
I haven't done so yet.  Let me know if this is a feature you would like.

\subsubsection{Input/Output}

The simplest output is the usual:
\begin{verbatim}
os << m
\end{verbatim}
where \tt{os} is any \tt{std::ostream}.
The output format is the same as for a \tt{Matrix}, including all the 0's.

There is also a compact format:
\begin{verbatim}
m.WriteCompact(os)
\end{verbatim}
which outputs in the format:
\begin{verbatim}
D n ( m(0,0)  m(1,1)  m(2,2)  ...  m(n-1,n-1) )
\end{verbatim}

The same (compact, that is) format can be read back in the usual two ways:
\begin{verbatim}
tmv::DiagMatrix<T> d(n);
is >> d;
std::auto_ptr<tmv::DiagMatrix<T> > pd;
is >> pd;
\end{verbatim}

One can write small values as 0 with
\begin{verbatim}
m.Write(std::ostream& os, RT thresh)
m.WriteCompact(std::ostream& os, RT thresh)
\end{verbatim}

\subsection{Upper/lower triangle matrices}

The \tt{UpperTriMatrix} class is our upper triangle matrix class, which is non-zero
only on the main diagonal and above.  \tt{LowerTriMatrix} is our class for lower
triangle matrices, which are non-zero only on the main diagonal and below.

The class \tt{UpperTriMatrix} inherits from \tt{GenUpperTriMatrix}, 
and the class \tt{LowerTriMatrix}
inherits from \tt{GenLowerTriMatrix}, both of
which in turn inherit from \tt{BaseMatrix}.  
The various views and composite classes described below
also inherit from \tt{GenUpperTriMatrix} and \tt{GenLowerTriMatrix} 
as appropriate.

All of the routines are analogous for \tt{UpperTriMatrix} and 
\tt{LowerTriMatrix}, so we only list each routine once 
(the \tt{UpperTriMatrix} version for
definiteness).  

The \tt{UpperTriMatrix} and \tt{LowerTriMatrix} routines are included by:
\begin{verbatim}
#include "TMV_Tri.h"
\end{verbatim}

In addition to the \tt{T} template parameter, there are three other template 
parameters:
\tt{dt} which can be either \tt{tmv::UnitDiag} or \tt{tmv::NonUnitDiag},
\tt{stor} which can be \tt{tmv::RowMajor} or \tt{tmv::ColMajor},
and \tt{index} which can be \tt{tmv::CStyle} or \tt{tmv::FortranStyle}.
The default values for these template parameters are 
\tt{NonUnitDiag}, \tt{RowMajor}, and \tt{CStyle} respectively.

The storage of both an \tt{UpperTriMatrix} and a \tt{LowerTriMatrix} takes
$N \times N$ elements of memory, even though approximately half of them 
are never used.  Someday, I'll write the packed storage versions which allow for
efficient storage of the matrices.

\subsubsection{Constructors}

\begin{itemize}
\item 
\begin{verbatim}
tmv::UpperTriMatrix<T,dt,stor,index> m(size_t n)
\end{verbatim}
Makes an \tt{n} $\times$ \tt{n} \tt{UpperTriMatrix} with \underline{uninitialized} values.
If debugging is turned on (i.e. not turned off
with -DNDEBUG), then the values are in fact initialized to 888.  

\item
\begin{verbatim}
tmv::UpperTriMatrix<T,dt,stor,index> m(size_t n, T x)
\end{verbatim}
Makes an \tt{n} $\times$ \tt{n} \tt{UpperTriMatrix} with all values equal to \tt{x}.

\item 
\begin{verbatim}
tmv::UpperTriMatrix<T,dt,stor,index> m(const GenMatrix<T>& mm)
tmv::UpperTriMatrix<T,dt,stor,index> m(const GenUperTriMatrix<T>& mm)
\end{verbatim}
Makes an \tt{UpperTriMatrix} which copies the corresponding values of
\tt{mm}.  Note that the second one is allowed to have \tt{mm} be
\tt{NonUnitDiag} but \tt{dt = UnitDiag}, in which case only the 
off-diagonal elements are copied.  The converse would set the 
diagonal of the new \tt{UpperTriMatrix} to all 1's.

\item
\begin{verbatim}
tmv::UpperTriMatrixView<T,index> m = 
        UpperTriMatrixViewOf(Matrix<T>& mm, DiagType dt)
tmv::UpperTriMatrixView<T,index> m = 
        UpperTriMatrixViewOf(const MatrixView<T>& mm, DiagType dt)
tmv::ConstUpperTriMatrixView<T,index> m = 
        UpperTriMatrixViewOf(const GenMatrix<T>& mm, DiagType dt)
tmv::UpperTriMatrixView<T,index> m = 
        UpperTriMatrixViewOf(UpperTriMatrix<T>& mm, DiagType dt)
tmv::UpperTriMatrixView<T,index> m = 
        UpperTriMatrixViewOf(const UpperTriMatrixView<T>& mm, 
        DiagType dt)
tmv::ConstUpperTriMatrixView<T,index> m = 
        UpperTriMatrixViewOf(const GenUpperTriMatrix<T>& mm, 
        DiagType dt)
\end{verbatim}
Make an \tt{UpperTriMatrixView} of the corresponding portion of \tt{mm}. 
Note the last three functions provide a way to review a 
\tt{NonUnitDiag} \tt{UpperTriMatrix} as \tt{UnitDiag}.

\item
\begin{verbatim}
tmv::UpperTriMatrixView<T,index> m = 
        tmv::UpperTriMatrixViewOf(T* mm, size_t n, DiagType dt, 
        StorageType stor)
tmv::ConstUpperTriMatrixView<T,index> m = 
        tmv::UpperTriMatrixViewOf(const T* mm, size_t n, DiagType dt, 
        StorageType stor)
\end{verbatim}
Make a \tt{UpperTriMatrixView} with the actual memory elements, 
\tt{mm} as the upper triangle.
One wrinkle here is that if \tt{dt} is \tt{UnitDiag}, then 
\tt{mm} is still the location of the
upper left corner, even though that value is never used 
(since the value is just taken to
be 1).  In fact, \tt{mm} must be of length \tt{n} $\times$ \tt{n},
so all of the lower triangle
elements must be in memory, even though they are never used.

\item
\begin{verbatim}
tmv::UpperTriMatrix<T> m1 = m2
tmv::UpperTriMatrix<T> m1(const GenUpperTriMatrix<T2>& m2)
\end{verbatim}
Copy the \tt{UpperTriMatrix m2}, which may be of any type \tt{T2} so long
as values of type \tt{T2} are convertible into type \tt{T}.
The assignment operator has the same flexibility.

\end{itemize}

\subsubsection{Access}

\begin{verbatim}
m.nrows() = m.ncols() = m.colsize() = m.rowsize() = m.size()
m(i,j)
m.row(size_t i, size_t j1, size_t j2)
m.col(size_t i, size_t j1, size_t j2)
m.diag()
m.diag(int i)
m.diag(int i, size_t k1, size_t k2)
\end{verbatim}
Note that the versions of \tt{row} and \tt{col} with only one argument are
missing, since the full row or column isn't accessible as a \tt{VectorView}.
You must specify a valid range within the row or column that you want, 
given the upper or lower triangle shape of \tt{m}.
If \tt{dt} is \tt{UnitDiag}, then the range may not include the diagonal element.
Similarly, \tt{m.diag()} is valid only if \tt{dt} is \tt{NonUnitDiag}.

\begin{verbatim}
m.SubVector(size_t i, size_t j, int istep, int jstep, size_t size)
m.SubMatrix(int i1, int i2, int j1, int j2)
m.SubMatrix(int i1, int i2, int j1, int j2, int istep, int jstep)
\end{verbatim}
This works the same as for \tt{Matrix}, except that all of the elements in the 
subvector or submatrix must be completely within the upper or lower triangle, as
appropriate.  If \tt{dt} is \tt{UnitDiag}, then no elements may be on the 
main diagonal.

\begin{verbatim}
m.SubTriMatrix(int i1, int i2, int istep = 1)
\end{verbatim}
This returns the upper or lower triangle matrix whose upper-left
corner is \tt{m(i1,i1)}, and whose lower-right corner is 
\tt{m(i2-istep,i2-istep)} for C-style indexing or \tt{m(i2,i2)} 
for Fortran-style indexing.  If \tt{istep} $\neq 1$, then it is the 
step in both the \tt{i} and \tt{j} directions.

\begin{verbatim}
m.OffDiag()
\end{verbatim}
This returns a view to the portion of the triangle matrix that does not include
the diagonal elements.  It will always be \tt{NonUnitDiag}.
Internally, it provides an easy way to deal with the \tt{UnitDiag} triangle matrices
for many routines.  But it may be useful for some users as well.

\begin{verbatim}
m.Real()
m.Imag()
m.Transpose()
m.Conjugate()
m.Adjoint()
m.View()
\end{verbatim}
Note that the transpose and adjoint of a \tt{LowerTriMatrix} is an
\tt{UpperTriMatrixView} and vice versa.

\subsubsection{Functions}

\begin{verbatim}
RT m.Norm1() = Norm1(m)
RT m.Norm2() = Norm2(m)
RT m.NormInf() = NormInf(m)
RT m.NormF() = NormF(m) = m.Norm() = Norm(m)
RT m.NormSq() = NormSq(m)
RT m.MaxAbsElement() = MaxAbsElement(m)
T m.Trace() = Trace(m)
T m.Det() = Det(m)
minv = m.Inverse() = Inverse(m)
m.Inverse(Matrix<T>& minv)
m.Inverse(UpperTriMatrix<T>& minv)
m.InverseATA(Matrix<T>& cov)
\end{verbatim}
Since the inverse of an \tt{Upper[Lower]TriMatrix} is also an 
\tt{Upper[Lower]TriMatrix},
we also provide a version of the \tt{Inverse} syntax which allows \tt{minv}
to be an \tt{Upper[Lower]TriMatrix}.  
The same option is available with the operator version (\tt{minv = m.Inverse()}).

\begin{verbatim}
m.Zero()
m.SetAllTo(T x)
m.Clip(RT thresh)
m.SetToIdentity(T x = 1)
m.ConjugateSelf()
Swap(m1,m2)
\end{verbatim}

\subsubsection{Arithmetic}

In addition to \tt{x}, \tt{v}, and \tt{m} from before, 
we now add \tt{U} and \tt{L} for a \tt{UpperTriMatrix}
and \tt{LowerTriMatrix} respectively.  Where the syntax is identical
for the two cases, only the \tt{U} form is listed.

\begin{verbatim}
U2 = -U1

U2 = x * U1
U2 = U1 [*/] x

U3 = U1 [+-] U2
m2 = m1 [+-] U
m2 = U [+-] m1
m = L [+-] U
m = U [+-] L

U [*/]= x

U2 [+-]= U1
m [+-]= U

v2 = U * v1
v2 = v1 * U
v *= U

U3 = U1 * U2
m2 = U * m1
m2 = m1 * U
m = U * L
m = L * U
U2 *= U1
m *= U

U2 = U1 [+-] x
U2 = x [+-] U1
U [+-]= x

\end{verbatim}

\subsubsection{Division}

The division operations are: (again omitting the L forms when redundant)
\begin{verbatim}
v2 = v1 [/%] U
m2 = m1 [/%] U
m2 = U [/%] m1
U3 = U1 [/%] U2
m = U [/%] L
m = L [/%] U
v [/%]= U
U2 [/%]= U1
m [/%]= U
\end{verbatim}

There is only one allowed \tt{DivType} for an \tt{UpperTriMatrix} or a 
\tt{LowerTriMatrix}: \tt{LU}.
And, since it is also the default behavior,
there is no reason to ever specify it.
Furthermore, as with a \tt{DiagMatrix},
the decomposition requires no work at all.  In fact, the ease of dividing by a 
upper or lower triangle matrix is precisely why the LU decomposition is useful.
Hence, it is always done in place.  i.e. no extra storage is needed, and all of
the \tt{m.DivideInPlace()}, \tt{m.SaveDiv()}, etc. are irrelevant.

If an \tt{UpperTriMatrix} or \tt{LowerTriMatrix} is singular, 
you can check easily with \tt{m.Singular()},
but there is no direct way to treat the division as a SVD and skip any
divisions by 0.  So trying to divide by a singular triangle matrix 
will result in a runtime error.

\subsubsection{Input/Output}

The simplest output is the usual:
\begin{verbatim}
os << m
\end{verbatim}
where \tt{os} is any \tt{std::ostream}.
The output format is the same as for a \tt{Matrix}, including all the 0's.

There is also a compact format.  For an \tt{UpperTriMatrix},
\begin{verbatim}
U.WriteCompact(os)
\end{verbatim}
outputs in the format:
\begin{verbatim}
U n 
( m(0,0)  m(0,1)  m(0,2)  ...  m(0,n-1) )
( m(1,1)  m(1,2)  ...  m(1,n-1) )
...
( m(n-1,n-1) )
\end{verbatim}
For a \tt{LowerTriMatrix},
\begin{verbatim}
L.WriteCompact(os)
\end{verbatim}
outputs in the format:
\begin{verbatim}
L n 
( m(0,0) )
( m(1,0)  m(1,1) )
...
( m(n-1,0)  m(n-1,1) ... m(n-1,n-1) )
\end{verbatim}

In each case, the compact format can be read back in the usual two ways:
\begin{verbatim}
tmv::UpperTriMatrix<T> U(n);
tmv::LowerTriMatrix<T> L(n);
is >> U >> L;
std::auto_ptr<tmv::UpperTriMatrix<T> > pU;
std::auto_ptr<tmv::LowerTriMatrix<T> > pL;
is >> pU >> pL;
\end{verbatim}

One can write small values as 0 with
\begin{verbatim}
m.Write(std::ostream& os, RT thresh)
m.WriteCompact(std::ostream& os, RT thresh)
\end{verbatim}

\subsection{Symmetric and hermitian matrices}

The \tt{SymMatrix} class is our symmetric matrix class.  A symmetric matrix is
one for which $m = m^T$.  We also have a class called
\tt{HermMatrix}, which is our hermitian matrix class.  A hermitian matrix
is one for which $m = m^\dagger$.  The two are exactly the same 
if \tt{T} is real, but for complex \tt{T}, they are different.

Both classes inherit from \tt{GenSymMatrix},
which has an internal parameter to keep track
of whether it is symmetric or hermitian.
\tt{GenSymMatrix} in turn inherits from \tt{BaseMatrix}.  
The various views and composite classes described below 
also inherit from \tt{GenSymMatrix}.

Usually, the symmetric/hermitian question does not affect the use of the classes.
(It does affect the actual calculations performed of course.)  So we will just refer
here to \tt{SymMatrix} and point out whenever a \tt{HermMatrix} acts differently.

One general caveat about complex \tt{HermMatrix} calculations is that the diagonal
elements should all be real.  Some calculations assume this, and only 
use the real part of the diagonal elements.  Other calculations use the 
full complex value as it is in memory.  Therefore, if you set a diagonal element 
to a non-real value at some point, the results will likely be wrong in
unpredictable ways.  Plus of course, your matrix won't actually be hermitian any more,
so the right answer is undefined in any case.

All the \tt{SymMatrix} and \tt{HermMatrix} routines are included by:
\begin{verbatim}
#include "TMV_Sym.h"
\end{verbatim}

In addition to the \tt{T} template parameter, there are three other template 
parameters: \tt{uplo} which can be either \tt{tmv::Upper} or \tt{tmv::Lower},
\tt{stor} which can be either \tt{tmv::RowMajor} or \tt{tmv::ColMajor}, and
\tt{index} which can be either \tt{tmv::CStyle} or \tt{tmv::FortranStyle}.
The parameter \tt{uplo} refers to which triangle the data are actually stored in, 
since the
other half of the values are identical, so we do not need to reference them.
The default values for these are \tt{Upper}, \tt{RowMajor}, and \tt{CStyle}
respectively.

The storage of a \tt{SymMatrix} takes
$N \times N$ elements of memory, even though approximately half of them 
are never used.  Someday, I'll write the packed storage versions which allow for
efficient storage of the matrices.

\subsubsection{Constructors}

\begin{itemize}
\item 
\begin{verbatim}
tmv::SymMatrix<T,uplo,stor,index> m(size_t n)
\end{verbatim}
Makes an \tt{n} $\times$ \tt{n} \tt{SymMatrix} with \underline{uninitialized} values.
If debugging is turned on (i.e. not turned off
with -DNDEBUG), then the values are in fact initialized to 888.  

\item
\begin{verbatim}
tmv::SymMatrix<T,uplo,stor,index> m(size_t n, T x)
\end{verbatim}
Makes an \tt{n} $\times$ \tt{n} \tt{SymMatrix} with all values equal to \tt{x}.
For the \tt{HermMatrix} version of this, \tt{x} must be real.

\item
\begin{verbatim}
tmv::SymMatrix<T,uplo,stor,index> m(size_t n, const T* mm)
tmv::SymMatrix<T,uplo,stor,index> m(size_t n, const std::vector<T>& mm)
\end{verbatim}
Makes a \tt{SymMatrix} which copies the elements from \tt{mm}.  

\item 
\begin{verbatim}
tmv::SymMatrix<T,uplo,stor,index> m(const GenMatrix<T>& mm)
\end{verbatim}
Makes a \tt{SymMatrix} which copies the corresponding values of \tt{mm}.  

\item
\begin{verbatim}
tmv::SymMatrixView<T,index> m = 
        tmv::SymMatrixViewOf(Matrix<T>& mm, UpLoType uplo)
tmv::SymMatrixView<T,index> m = 
        tmv::SymMatrixViewOf(const MatrixView<T>& mm,UpLoType uplo)
tmv::ConstSymMatrixView<T,index> m = 
        tmv::SymMatrixViewOf(const GenMatrix<T>& mm, UpLoType uplo)
tmv::SymMatrixView<T,index> m = 
        tmv::HermMatrixViewOf(Matrix<T>& mm, UpLoType uplo)
tmv::SymMatrixView<T,index> m = 
        tmv::HermMatrixViewOf(const MatrixView<T>& mm, UpLoType uplo)
tmv::ConstSymMatrixView<T,index> m = 
        tmv::HermMatrixViewOf(const GenMatrix<T>& mm, UpLoType uplo)
\end{verbatim}
Make a \tt{SymMatrixView} of the corresponding portion of \tt{mm}.  

\item
\begin{verbatim}
tmv::SymMatrixView<T,index> m = 
        tmv::SymMatrixViewOf(T* mm, size_t n, UpLoType uplo, 
        StorageType stor)
tmv::ConstSymMatrixView<T,index> m = 
        tmv::SymMatrixViewOf(const T* mm, size_t n, UpLoType uplo, 
        StorageType stor)
tmv::SymMatrixView<T,index> m = 
        tmv::HermMatrixViewOf(T* mm, size_t n, UpLoType uplo, 
        StorageType stor)
tmv::ConstSymMatrixView<T,index> m = 
        tmv::HermMatrixViewOf(const T* mm, size_t n, UpLoType uplo, 
        StorageType stor)
\end{verbatim}
Make a \tt{SymMatrixView} with the actual memory elements, \tt{mm} in either the 
upper or lower triangle.
\tt{mm} must be of length \tt{n} $\times$ \tt{n}, even though only about half 
of the values are actually used,

\item
\begin{verbatim}
tmv::SymMatrix<T,uplo,stor,index> m1 = m2
tmv::SymMatrix<T,uplo,stor,index> m1(
        const GenSymMatrix<T2,uplo2,stor2,index2>& m2)
\end{verbatim}
Copy the \tt{SymMatrix m2}, which may be of any type \tt{T2} so long
as values of type \tt{T2} are convertible into type \tt{T}.
The assignment operator has the same flexibility.
If \tt{T} and \tt{T2} are complex, then \tt{m1} and \tt{m2} need to be
either both symmetric or both hermitian.

\end{itemize}

\subsubsection{Access}

\begin{verbatim}
m.nrows() = m.ncols() = m.colsize() = m.rowsize() = m.size()
m(i,j)
m.row(size_t i, size_t j1, size_t j2)
m.col(size_t i, size_t j1, size_t j2)
m.diag()
m.diag(int i)
m.diag(int i, size_t k1, size_t k2)
\end{verbatim}
As for triangle matrices, the versions of \tt{row} and \tt{col} with only one argument are
missing, since the full row or column isn't accessible as a \tt{VectorView}.
You must specify a valid range within the row or column that you want, 
given the upper or lower storage of \tt{m}.
The diagonal element may be in a \tt{VectorView} with either elements in the 
lower triangle or the upper triangle, but not both.  To access a full row, you would 
therefore need to use two steps:
\begin{verbatim}
m.row(i,0,i) = ...
m.row(i,i,ncols) = ...
\end{verbatim}

\begin{verbatim}
m.SubVector(size_t i, size_t j, int istep, int jstep, size_t size)
m.SubMatrix(int i1, int i2, int j1, int j2)
m.SubMatrix(int i1, int i2, int j1, int j2, int istep, int jstep)
\end{verbatim}
These work the same as for a \tt{Matrix}, except that the entire
subvector or submatrix must be completely within the upper or lower triangle.

\begin{verbatim}
m.SubSymMatrix(int i1, int i2, int istep = 1)
\end{verbatim}
This returns a \tt{SymMatrixView} of \tt{m} whose upper-left
corner is \tt{m(i1,i1)}, and whose lower-right corner is 
\tt{m(i2-istep,i2-istep)}.  If \tt{istep} $\neq 1$, then it is the 
step in both the \tt{i} and \tt{j} directions.

\begin{verbatim}
m.UpperTri(DiagType dt)
m.LowerTri(DiagType dt)
\end{verbatim}
These return an \tt{UpperTriMatrixView} or \tt{LowerTriMatrixView} respectively of the 
corresponding portion of the \tt{SymMatrix}.  Both of these are valid calls regardless
of which triangle stores the actual data for \tt{m}.  The argument \tt{dt} should
be either \tt{UnitDiag} or \tt{NonUnitDiag}, or it may be omitted, in which case
\tt{NonUnitDiag} is assumed.

\begin{verbatim}
m.Real()
m.Imag()
m.Transpose()
m.Conjugate()
m.Adjoint()
m.View()
\end{verbatim}
Note that the imaginary part of a complex hermitian matrix is anti-symmetric,
which is not yet supported as a special matrix type,
so \tt{m.Imag()} is illegal for a \tt{HermMatrix}.  If you need to
manipulate the imaginary part of a \tt{HermMatrix},
you could use \tt{m.UpperTri().Imag()}.  Or, 
since the diagonal elements are all real.
you could also use \tt{m.UpperTri().OffDiag().Imag()}. 

\subsubsection{Functions}

\begin{verbatim}
RT m.Norm1() = Norm1(m)
RT m.Norm2() = Norm2(m)
RT m.NormInf() = NormInf(m)
RT m.NormF() = NormF(m) = m.Norm() = Norm(m)
RT m.NormSq() = NormSq(m)
RT m.MaxAbsElement() = MaxAbsElement(m)
T m.Trace() = Trace(m)
T m.Det() = Det(m)
minv = m.Inverse() = Inverse(m)
m.Inverse(Matrix<T>& minv)
m.Inverse(SymMatrix<T>& minv)
m.InverseATA(Matrix<T>& cov)
\end{verbatim}
Since the inverse of a \tt{SymMatrix} is also a \tt{SymMatrix},
we also provide a version of the \tt{Inverse} syntax which allow minv
to be a \tt{SymMatrix}.  The same option is available with the operator version
(\tt{minv = m.Inverse()}).

\begin{verbatim}
m.Zero()
m.SetAllTo(T x)
m.Clip(RT thresh)
m.SetToIdentity(T x = 1)
m.ConjugateSelf()
m.TransposeSelf()
m.SwapRowsCols(size_t i1, size_t i2)
Swap(m1,m2)
\end{verbatim}
\tt{TransposeSelf} does nothing to a \tt{SymMatrix}, and it is equivalent to
\tt{ConjugateSelf} for a \tt{HermMatrix}.
The new method, \tt{SwapRowsCols}, would be equivalent to 
\begin{verbatim}
m.SwapRows(i1,i2).SwapCols(i1,i2);
\end{verbatim}
except that neither of these functions are allowed for a \tt{SymMatrix}, since 
they result in non-symmetric matrices.  Only the combination of both
maintains the symmetry of the matrix.  So this combination is included as
a method.

\subsubsection{Arithmetic}
\label{symarith}

In addition to \tt{x}, \tt{v}, and \tt{m} from before,
we now add \tt{s} for a \tt{SymMatrix}.

\begin{verbatim}
s2 = -s1

s2 = x * s1
s2 = s1 [*/] x

s3 = s1 [+-] s2
m2 = m1 [+-] s
m2 = s [+-] m1

s [*/]= x

s2 [+-]= s1
m [+-]= s

v2 = s * v1
v2 = v1 * s
v *= s

m = s1 * s2
m2 = s * m1
m2 = m1 * s
m *= s

s2 = s1 [+-] x
s2 = x [+-] s1
s [+-]= x

s = v ^ v
s [+-]= v ^ v
s = m * m.Transpose()
s [+-]= m * m.Transpose()
s = U * U.Transpose()
s [+-]= U * U.Transpose()
s = L * L.Transpose()
s [+-]= L * L.Transpose()
\end{verbatim}
For outer products, both \tt{v}'s need to be the same actual data.  If \tt{s}
is complex hermitian, then it should actually be \tt{s = v \^{} v.Conjugate}.
Likewise for the next one (called a ``rank-k update''), the \tt{m}'s need to be the
same data, and for a complex hermitian \tt{s}, \tt{Transpose}
should be replaced with \tt{Adjoint}.

There is a minor issue with \tt{HermMatrix} arithmetic that may occasionally be
worth knowing about.
Some operations with a \tt{HermMatrix}, as far as the compiler can tell, may or may
not result in another \tt{HermMatrix}.  For example, a scalar times a \tt{HermMatrix}
is hermitian if the scalar is real, but not if the scalar is complex.
Likewise for the sum of a \tt{HermMatrix} and a scalar.
Also, since my base class is the same for \tt{SymMatrix} and \tt{HermMatrix},
the sum of two \tt{GenSymMatrix}es may or may not be symmetric or 
hermitian.

So the \tt{SymMatrixComposite} class, which is the return type of all of these operations are derived from \tt{GenMatrix},
rather than from \tt{GenSymMatrix}.  This means that if you let it
self-instantiate, rather than assign it to a \tt{SymMatrix}, it will 
instantiate as a regular \tt{Matrix}.  Most of the time, this won't matter, since
you will generally assign the result to either a \tt{SymMatrix} or a \tt{Matrix}
as appropriate.  However, some things that should be allowed are not for this
reason.  For example:
\begin{verbatim}
s3 += x*s1+s2;
\end{verbatim}
is not legal, even if everything is real.  This is because there is no 3-matrix 
\tt{Add} function.  (The \tt{AddMM} function below just adds a multiple of one 
matrix to another.)  So the right hand side needs to be instantiated before
being added to the left side, and it will instantiate as a regular \tt{Matrix},
which cannot be added to a \tt{SymMatrix}.  

One work-around is to explicitly tell the compiler to instantiate the right hand
side as a \tt{SymMatrix}:
\begin{verbatim}
s3 += SymMatrix<T>(x*s1+s2);
\end{verbatim}
If the ``\tt{+=}'' had been just ``\tt{=}'', then we wouldn't have had to do this,
since the composite object which stores \tt{x*s1+s2} is assignable to a \tt{SymMatrix}
despite being derived from \tt{GenMatrix}.

Another work-around, which I suspect will usually be preferred, is to break the 
equation into multiple statements, each of which are simple enough to not require
any instantiation:
\begin{verbatim}
s3 += x*s1;
s3 += s2;
\end{verbatim}

\subsubsection{Division}

The division operations are:
\begin{verbatim}
v2 = v1 [/%] s
m2 = m1 [/%] s
m2 = s [/%] m1
m = s1 [/%] s2
v [/%]= s
m [/%]= s
\end{verbatim}

\tt{SymMatrix} has three possible choices for the division decomposition:
\begin{enumerate}
\item
\tt{m.DivideUsing(tmv::LU)} will perform something similar to the LU decomposition
for regular matrices.  Instead, it does what is called an LDL or Bunch-Kauffman
decomposition.  

A permutation of \tt{m} is decomposed into a lower triangle matrix ($L$)
times a symmetric block diagonal matrix ($D$) times the transpose of $L$.
$D$ has either 1x1 and 2x2 blocks down the diagonal.  For hermitian matrices,
the third term is the adjoint of $L$ rather than the transpose.

This is the default decomposition to use if you don't specify anything.

To access this decomposition, use:
\begin{verbatim}
ConstLowerTriMatrixView<T> m.LUD().GetL()
Matrix<T> m.LUD().GetD()
size_t* m.LUD().GetP()
\end{verbatim}
The following should result in a matrix numerically very close to \tt{m}.
\begin{verbatim}
Matrix<T> m2 = m.LUD().GetL() * m.LUD().GetD() * 
        m.LUD().GetL().Transpose();
m2.ReversePermuteRows(m.LUD().GetP());
m2.ReversePermuteCols(m.LUD().GetP());
\end{verbatim}
For a complex, hermitian \tt{m}, you would need to replace
\tt{Transpose} with \tt{Adjoint}.

\item
\tt{m.DivideUsing(tmv::CH)} will perform a Cholesky decomposition.  
\tt{m} must be hermitian (or real symmetric) to use \tt{CH}, since that is the
only kind of matrix that has a Cholesky decomposition.  

It is also similar to an 
LU decomposition, where $U$ is the adjoint of $L$, and there is no permutation.
It can be a bit dangerous, since not all hermitian matrices have such a decomposition,
so the decomposition could fail.  Only so-called ``positive-definite'' hermitian 
matrices have a Cholesky decomposition.  A positive-definite matrix has
all positive real eigenvalues.  In general, hermitian matrices have real, but
not necessarily positive eigenvalues.  

One example of a 
positive-definite matrix is $m = A^\dagger A$ where $A$ is any matrix.   
Then $m$ is guaranteed to be positive-semi-definite
(which means some of the eigenvalues may be 0, but not negative).
In this case, the routine will usually work, but still might fail from 
numerical round-off errors if $m$ is nearly singular.  
Even if it doesn't fail, if \tt{m} is nearly singular, the
calculation may be numerically unstable, and SVD may be a better choice.

The only advantage of Cholesky over Bunch-Kauffman is speed.  If you know your 
matrix is positive-definite, the Cholesky decomposition is the fastest way to 
do division.

To access this decomposition, use:
\begin{verbatim}
ConstLowerTriMatrixView<T> m.CHD().GetL()
\end{verbatim}
The following should result in a matrix numerically very close to \tt{m}.
\begin{verbatim}
Matrix<T> m2 = m.CHD().GetL() * m.CHD().GetL().Adjoint()
\end{verbatim}

\item
\tt{m.DivideUsing(tmv::SV)} will perform a singular value decomposition.

For hermitian matrices (including real symmetric matrices), 
the singular value decomposition is quite elegant,
since it is also essentially an eigenvalue decomposition.
In this case $V = U^\dagger$, $S$ is a diagonal matrix of the eigenvalues,
and the columns of $U$ are the corresponding eigenvectors.

The one slight change in how we make the decomposition
is that we have to allow the elements of $S$ to be negative, since eigenvalues of 
a matrix can be negative, whereas singular 
values are really defined to be positive.  So if you want to use them as
singular values, you just need to remember to take the absolute value.
The eigenvalue decomposition has important uses in 
its own right beyond just division, so we didn't want to spoil that.

To access this decomposition, use:
\begin{verbatim}
ConstMatrixView<T> m.SVD().GetU()
ConstVectorView<RT> m.SVD().GetS()
ConstMatrixView<T> m.SVD().GetV()
\end{verbatim}
The following should result in a matrix numerically very close to \tt{m}.
\begin{verbatim}
Matrix<T> m2 = m.SVD().GetU() * DiagMatrixViewOf(m.SVD().GetS()) * 
        m.SVD().GetV()
\end{verbatim}

For a complex, symmetric \tt{m}, the situation is not nearly so elegant.
Technically, one could find a decomposition $m = USU^T$ where $S$ is 
a complex diagonal matrix, but such a decomposition is not easy to find. 
Standard SVD algorithms find $S$ with real values.  In this case 
$V$ is not $U^T$.  So the algorithm for complex symmetric matrices 
just finds $m = USV$, although it does use the symmetry of the matrix to 
speed up portions of the algorithm relative that that for a generic matrix.

The access is also necessarily different, since the object returned by 
\tt{m.SVD()} implicitly assumes that \tt{V = U.Adjoint()}, so we need a 
new accessor: \tt{m.SymSVD()}.  Other than that, the access syntax is 
the same as for hermitian matrices.

Both versions also have the same control and access routines as a regular SVD:
\begin{verbatim}
m.SVD().Thresh(RT thresh)
m.SVD().Top(size_t nsing)
RT m.SVD().Condition()
size_t m.SVD().GetKMax()
\end{verbatim}
(Likewise for m.SymSVD().)

There are also the \tt{SVS}, \tt{SVU}, and \tt{SVV} options for the decomposition
which cannot be used for division, but which can access parts of the decomposition.

\end{enumerate}

For all three \tt{DivType}s, the routines 
\begin{verbatim}
m.SaveDiv()
m.SetDiv()
m.ReSetDiv()
m.UnSetDiv()
m.DivideInPlace()
\end{verbatim}
work the same as for regular \tt{Matrix}es.

\subsubsection{Input/Output}

The simplest output is the usual:
\begin{verbatim}
os << m
\end{verbatim}
where \tt{os} is any \tt{std::ostream}.
The output format is the same as for a \tt{Matrix}.

There is also a compact format for a \tt{SymMatrix}:
\begin{verbatim}
m.WriteCompact(os)
\end{verbatim}
outputs in the format:
\begin{verbatim}
H/S n 
( m(0,0) )
( m(1,0)  m(1,1) )
...
( m(n-1,0)  m(n-1,1) ... m(n-1,n-1) )
\end{verbatim}
where \tt{H/S} means \underline{either} the character \tt{H} or \tt{S}, which indicates whether 
the matrix is hermitian or symmetric.

In each case, the same compact format can be read back in the usual two ways:
\begin{verbatim}
tmv::SymMatrix<T> s(n);
tmv::HermMatrix<T> h(n);
is >> s >> h;
std::auto_ptr<tmv::SymMatrix<T> > ps;
std::auto_ptr<tmv::HermMatrix<T> > ph;
is >> ps >> ph;
\end{verbatim}

One can write small values as 0 with
\begin{verbatim}
m.Write(std::ostream& os, RT thresh)
m.WriteCompact(std::ostream& os, RT thresh)
\end{verbatim}

\subsection{Band matrices}

The \tt{BandMatrix} class is our band-diagonal matrix, which is only non-zero
on the main diagonal and a few sub- and super-diagonals.  While band-diagonal
matrices are usually square, we allow for non-square banded matrices as well.
You may even have rows or columns which are completely outside of the 
band structure, and hence are all 0.  For example a $10\times 5$ band matrix
with 2 sub-diagonals is valid even though the bottom 3 rows are all 0.

Throughout, we use \tt{nlo} to refer to the number of sub-diagonals 
(below the main diagonal) stored
in the \tt{BandMatrix}, and \tt{nhi} to refer to the number of super-diagonals
(above the main diagonal).

\tt{BandMatrix} inherits from \tt{GenBandMatrix}, which in turn inherits from
\tt{BaseMatrix}.  
The various views and composite classes described below 
also inherit from \tt{GenBandMatrix}.

All the \tt{BandMatrix} routines are included by:
\begin{verbatim}
#include "TMV_Band.h"
\end{verbatim}

In addition to the \tt{T} template parameter, we also have \tt{stor} to indicate 
which storage you want to use, along with the usual \tt{index} parameter.
For this class, we add an additional storage possibility - 
along with \tt{RowMajor} and 
\tt{ColMajor} storage, a \tt{BandMatrix} may also be \tt{DiagMajor}, which 
has unit step along the diagonals.
The default values for these are \tt{RowMajor} and \tt{CStyle}.

For each type of storage, we require that the step size in each direction
be uniform within a given row, column or diagonal. 
This means that we require a few extra elements
of memory which are not actually used.

To demonstrate the different storage orders and why extra memory is required, 
here are three $6 \times 6$ band-diagonal
matrices, each with \tt{nlo = 2} and \tt{nhi = 3} in each of the different storage
types.  The number in each place indicates the offset in memory from the 
top left element.

\begin{align*}
\textrm{ColMajor:} ~ ~ & \left(\begin{array}{cccccc}0 & 5 & 10 & 15 &  &  \\1 & 6 & 11 & 16 & 21 &  \\2 & 7 & 12 & 17 & 22 & 27 \\ & 8 & 13 & 18 & 23 & 28 \\ &  & 14 & 19 & 24 & 29 \\ &  &  & 20 & 25 & 30\end{array}\right) \\
\textrm{RowMajor:} ~  ~ & \left(\begin{array}{cccccc}0 & 1 & 2 & 3 &  &  \\ 5 & 6 & 7 & 8 & 9 &  \\ 10 & 11 & 12 & 13 & 14 & 15 \\ & 16 & 17 & 18 & 19 & 20 \\ &  & 22 & 23 & 24 & 25 \\ &  &  & 28 & 29 & 30\end{array}\right) \\
\textrm{DiagMajor:} ~  ~ & \left(\begin{array}{cccccc}0 & 6 & 12 & 18 &  &  \\-5 & 1 & 7 & 13 & 19 &  \\ -10 & -4 & 2 & 8 & 14 & 20 \\ & -9 & -3 & 3 & 9 & 15 \\ &  & -8 & -2 & 4 & 10 \\ &  &  & -7 & -1 & 5\end{array}\right)
\end{align*}

First, notice that all three storage methods require 4 extra locations in memory 
which do not hold
any actual matrix data.  
(They require a total of 31 memory addresses for only 27 that are used.)
This is because we want to have the same step size between consecutive row elements
for every row.  Likewise for the columns (which in turn implies that it is also 
true for the diagonals).

For $N\times N$ square matrices, the total memory needed is $(N-1)*(nlo+nhi+1)+1$, 
which wastes $(nlo-1)*nlo/2 + (nhi-1)*nhi/2$ locations.  For non-square
matrices, the formula is more complicated, and changes slightly between the 
three storages.
If you want to know the memory used by a \tt{BandMatrix}, we provide the routine:
\begin{verbatim}
size_t BandStorageLength(StorageType stor, size_t nrows, size_t ncols, 
        int nlo, int nhi)
\end{verbatim}
  
For square matrices, all three methods always need the same amount of memory 
(and for non-square, they aren't very different), so the 
decision about which method to use should generally be based on performance 
considerations rather than memory usage.
The speed of the various matrix operations are different for the different storage 
types.  If the matrix calculation speed is important, it may be worth trying 
all three to see which is fastest for the operations you are using.

Second, notice that the \tt{DiagMajor} storage doesn't start with the 
upper left element as usual.
Rather, it starts at the start of the lowest sub-diagonal.  
So for the constructors which 
take the matrix information from an array (\tt{T*} or \tt{vector<T>}), 
the start of the array needs to be at the start of the lowest sub-diagonal.  

\subsubsection{Constructors}

\begin{itemize}
\item 
\begin{verbatim}
tmv::BandMatrix<T,stor,index> m(size_t nrows, size_t ncols, 
        int nlo, int nhi)
\end{verbatim}
Makes a \tt{BandMatrix} with \tt{nrows} rows, \tt{ncols} columns, 
\tt{nlo} sub-diagonals,
and \tt{nhi} super-diagonals with \underline{uninitialized} values.
If debugging is turned on (i.e. not turned off 
with -DNDEBUG), then the values are initialized to 888.

\item
\begin{verbatim}
tmv::BandMatrix<T,stor,index> m(size_t nrows, size_t ncols, 
        int nlo, int nhi, T x)
\end{verbatim}
Makes a \tt{BandMatrix} with all values equal to \tt{x}.

\item
\begin{verbatim}
tmv::BandMatrix<T,stor,index> m(size_t nrows, size_t ncols, 
        int nlo, int nhi, const T* mm)
tmv::BandMatrix<T,stor,index> m(size_t nrows, size_t ncols, 
        int nlo, int nhi, const std::vector<T>& mm)
\end{verbatim}
Makes a \tt{BandMatrix} which copies the elements from \tt{mm}.  See the 
discussion above about the different storage types to see what order these elements
should be.  The function \tt{BandStorageLength} will tell you how
long \tt{mm} must be.  The elements which don't fall in the bounds of the 
actual matrix are not used and may be left undefined.

\item 
\begin{verbatim}
tmv::BandMatrix<T,stor,index> m(const GenMatrix<T>& mm, 
        int nlo, int nhi)
tmv::BandMatrix<T,stor,index> m(const GenBandMatrix<T>& mm, 
        int nlo, int nhi)
\end{verbatim}
Makes an \tt{BandMatrix} which copies the corresponding values of \tt{mm}.  
For the second one, \tt{nlo} and \tt{nhi} must not be larger than those for \tt{mm}.

\item
\begin{verbatim}
tmv::BandMatrixView<T,index> m = 
        BandMatrixViewOf(Matrix<T>& mm, int nlo, int nhi)
tmv::BandMatrixView<T,index> m = 
        BandMatrixViewOf(const MatrixView<T>& mm, int nlo, int nhi)
tmv::ConstBandMatrixView<T,index> m = 
        BandMatrixViewOf(const GenMatrix<T>& m, int nlo, int nhi)
tmv::BandMatrixView<T,index> m = 
        BandMatrixViewOf(BandMatrix<T>& mm, int nlo, int nhi)
tmv::BandMatrixView<T,index> m = 
        BandMatrixViewOf(const BandMatrixView<T>& mm, int nlo, int nhi)
tmv::ConstBandMatrixView<T,index> m = 
        BandMatrixViewOf(const GenBandMatrix<T>& m, int nlo, int nhi)
\end{verbatim}
Make an \tt{BandMatrixView} of the corresponding portion of \tt{mm}.  

Note: if you want a \tt{BandMatrix} which includes the whole upper triangle
plus only a few sub-diagonals (or vice versa), 
you are better off memory-wise making a regular 
\tt{Matrix} and then viewing the portion you want as a \tt{BandMatrixView}.
This is because the memory requirements (for a square matrix anyway)
are identical for a regular \tt{Matrix} and a \tt{BandMatrix} when
$\tt{nlo} + \tt{nhi} = N$.  For a Hessenberg matrix (upper triangle plus 
one sub-diagonal),
this is already the case.  So with any more sub-diagonals, you are actually using
more memory with a \tt{BandMatrix} than with a \tt{Matrix}.  But you can still get the 
algorithmic speed-up of the \tt{BandMatrix} by viewing the portion you want
as a \tt{BandMatrixView} using the above constructor.

\item
\begin{verbatim}
tmv::BandMatrixView<T> m = 
        tmv::BandMatrixViewOf(T* mm, size_t ncols, size_t nrows, 
          int nlo, int nhi, StorageType stor)
tmv::ConstBandMatrixView<T> m = 
        tmv::BandMatrixViewOf(const T* mm, size_t ncols, size_t nrows, 
          int nlo, int nhi, StorageType stor)
\end{verbatim}
Make a \tt{BandMatrixView} using the actual memory elements, \tt{mm}.

\item
\begin{verbatim}
tmv::BandMatrix<T,stor> m = UpperBiDiagMatrix(const GenVector<T>& v1, 
        const GenVector<T>& v2)
tmv::BandMatrix<T,stor> m = LowerBiDiagMatrix(const GenVector<T>& v1, 
        const GenVector<T>& v2)
tmv::BandMatrix<T,stor> m = TriDiagMatrix(const GenVector<T>& v1,
        const GenVector<T>& v2, const GenVector<T>& v3)
\end{verbatim}
Shorthand to create bi- or tri-diagonal \tt{BandMatrix}es if you already have the 
\tt{Vector}s.  The \tt{Vector}s are in order from bottom to top in each case.

\item
\begin{verbatim}
tmv::BandMatrix<T> m1 = m2
tmv::BandMatrix<T> m1(const GenBandMatrix<T2>& m2)
\end{verbatim}
Copy the \tt{BandMatrix m2}, which may be of any type \tt{T2} so long
as values of type \tt{T2} are convertible into type \tt{T}.
The assignment operator has the same flexibility.

\end{itemize}

\subsubsection{Access}

\begin{verbatim}
m.nrows() = m.colsize()
m.ncols() = m.rowsize()
m(i,j)
m.row(size_t i, size_t j1, size_t j2)
m.col(size_t i, size_t j1, size_t j2)
m.diag()
m.diag(int i)
m.diag(int i, size_t k1, size_t k2)
\end{verbatim}
Again, the versions of \tt{row} and \tt{col} with only one argument are
missing, since the full row or column isn't accessible as a \tt{VectorView}.
You must specify a valid range within the row or column that you want, 
given the banded storage of \tt{m}.

\begin{verbatim}
m.SubVector(size_t i, size_t j, int istep, int jstep, size_t size)
m.SubMatrix(int i1, int i2, int j1, int j2)
m.SubMatrix(int i1, int i2, int j1, int j2, int istep, int jstep)
\end{verbatim}
These work the same as for a \tt{Matrix}, except that the entire
subvector or submatrix must be completely within the band.

\begin{verbatim}
m.SubBandMatrix(int i1, int i2, int j1, int j2, int newnlo, int newnhi)
m.SubBandMatrix(int i1, int i2, int j1, int j2, int newnlo, int newnhi, 
        int istep, int jstep)
\end{verbatim}
This returns a \tt{BandMatrixView} of a subset of a \tt{BandMatrix}.
The \tt{newnlo} and \tt{newnhi} parameters do not need to be different
from the existing \tt{nlo} and \tt{nhi} if \tt{i1} = \tt{j1} 
(i.e. the new main diagonal
is part of the old main diagonal).  However, if you are moving the upper left
corner off of the diagonal, you need to adjust \tt{nlo} and \tt{nhi} appropriately.
For example, if \tt{m} is a $6 \times 6$ \tt{BandMatrix} with 2 sub-diagonals and 
3 super-diagonals
(like our example above), the 3 super-diagonals may be viewed as a \tt{BandMatrixView}
with \tt{m.SubBandMatrix(0,5,1,6,0,2)}.

\begin{verbatim}
m.Rows(int i1,int i2)
m.Cols(int j1,int j2)
m.Diags(int k1, int k2)
\end{verbatim}
These return a \tt{BandMatrixView} of the parts of these rows,
columns or diagonals that
appear within the original banded structure.  For our example of viewing just
the super-diagonals of a $6 \times 6$ \tt{BandMatrix} with 2 sub- and 
3 super-diagonals, we
could instead use \tt{m.Diags(1,4)}.  The last 3 rows would be \tt{m.Rows(3,6)}.
Note that this wold be a $3 \times 5$ matrix with 0 sub-diagonals and 
4 super-diagonals.
These routines calculate the appropriate changes in the size and shape to include
all of the relevant parts of the rows or columns.

\begin{verbatim}
m.UpperBand()
m.LowerBand()
\end{verbatim}
These return a \tt{BandMatrixView} including the main diagonal and either the
super- or sub-diagonals, respectively.  As with the above methods, the size
is automatically set appropriately to include the entire band.  (This is only
non-trivial for non-square band matrices.)

\begin{verbatim}
m.Real()
m.Imag()
m.Transpose()
m.Conjugate()
m.Adjoint()
m.View()
\end{verbatim}

\subsubsection{Functions}

\begin{verbatim}
RT m.Norm1() = Norm1(m)
RT m.Norm2() = Norm2(m)
RT m.NormInf() = NormInf(m)
RT m.NormF() = NormF(m) = m.Norm() = Norm(m)
RT m.NormSq() = NormSq(m)
RT m.MaxAbsElement() = MaxAbsElement(m)
T m.Trace() = Trace(m)
T m.Det() = Det(m)
minv = m.Inverse() = Inverse(m)
m.Inverse(Matrix<T>& minv)
m.InverseATA(Matrix<T>& cov)
\end{verbatim}
The inverse of a \tt{BandMatrix} is not (in general) banded.  So \tt{minv} here
must be a regular \tt{Matrix}.

\begin{verbatim}
m.Zero()
m.SetAllTo(T x)
m.Clip(RT thresh)
m.SetToIdentity(T x = 1)
m.ConjugateSelf()
m.TransposeSelf()
Swap(m1,m2)
\end{verbatim}

\subsubsection{Arithmetic}

In addition to \tt{x}, \tt{v}, and \tt{m} from before, we now add \tt{b} 
for a \tt{BandMatrix}.

\begin{verbatim}
b2 = -b1

b2 = x * b1
b2 = b1 [*/] x

b3 = b1 [+-] b2
m2 = m1 [+-] b
m2 = b [+-] m1

b [*/]= x

b2 [+-]= b1
m [+-]= b

v2 = b * v1
v2 = v1 * b
v *= b

b3 = b1 * b2
m2 = b * m1
m2 = m1 * b
m *= b

b2 = b1 [+-] x
b2 = x [+-] b1
b [+-]= x
\end{verbatim}

\subsubsection{Division}

The division operations are:
\begin{verbatim}
v2 = v1 [/%] b
m2 = m1 [/%] b
m2 = b [/%] m1
m = b1 [/%] b2
v [/%]= b
m [/%]= b
\end{verbatim}

\tt{BandMatrix} has three possible choices for the division decomposition:
\begin{enumerate}
\item
\tt{m.DivideUsing(tmv::LU)} does a normal LU decomposition,
taking into account the band structure of the matrix which greatly speeds up 
the calculation into the lower and upper (banded) triangles.
This is the default decomposition to use for a square \tt{BandMatrix} 
if you don't specify anything.

This decomposition can only really
be done in place if either \tt{nlo} or \tt{nhi} is 0,
in which case it is automatically done in place,
since the \tt{BandMatrix} is already lower or upper triangle.
Thus, there is usually no reason to use the \tt{DivideInPlace()} method.

If this is not the case, and you really want to do the decomposition in place, you can
declare a matrix with a wider band and view the portion that represents the matrix
you actually want.  This view then can be divided in place.  More specifically,
you need to declare the wider \tt{BandMatrix} with \tt{ColMajor} storage,
with the smaller of
\{\tt{nlo}, \tt{nhi}\} as the number of sub-diagonals, and with 
(\tt{nlo} + \tt{nhi}) as the number of super-diagonals.  
Then you can use \tt{BandMatrixViewOf} to view the portion you want, transposing it if
necessary.  On the other hand, you are probably not going to get much of a speed gain 
from all of this finagling, so unless you are really memory starved, it's probably not worth it.

To access this decomposition, use:
\begin{verbatim}
bool m.LUD().IsTrans()
tmv::LowerTriMatrix<T,UnitDiag> m.LUD().GetL()
tmv::ConstBandMatrixView<T> m.LUD().GetU()
size_t* m.LUD().GetP()
\end{verbatim}
The following should result in a matrix numerically very close to \tt{m}.
\begin{verbatim}
tmv::Matrix<T> m2 = m.LUD().GetL() * m.LUD().GetU();
m2.ReversePermuteRows(m.LUD().GetP());
if (m.LUD().IsTrans()) m2.TransposeSelf();
\end{verbatim}

\item
\tt{m.DivideUsing(tmv::QR)} will perform a QR decomposition.  
This is the default method for a non-square \tt{BandMatrix}.

The same kind of convolutions need to be done to perform this in place as 
for the LU decomposition.

To access this decomposition, use:
\begin{verbatim}
bool m.QRD().IsTrans()
tmv::Matrix<T> m.QRD().GetQ()
tmv::ConstBandMatrixView<T> m.QRD().GetR()
\end{verbatim}
The following should result in a matrix numerically very close to \tt{m}.
\begin{verbatim}
tmv::Matrix<T> m2(m.nrows,m.ncols);
tmv::MatrixView<T> m2v = 
        m.QRD().IsTrans() ? m2.Transpose() : m2.View();
m2v = m.QRD().GetL() * m.QRD().GetU();
m2v.ReversePermuteRows(m.QRD().GetP());
\end{verbatim}

\item
\tt{m.DivideUsing(tmv::SV)} will perform a singular value decomposition.

This cannot be done in place.

To access this decomposition, use:
\begin{verbatim}
tmv::ConstMatrixView<T> m.BandSVD().GetU()
tmv::ConstVectorView<RT> m.BandSVD().GetS()
tmv::ConstMatrixView<T> m.BandSVD().GetV()
\end{verbatim}
The product of these three (using \tt{DiagMatrixViewOf(m.BandSVD().GetS())} 
for $S$)
should result in a matrix numerically very close to \tt{m}.

There are the same control and access routines as for a regular SVD:
\begin{verbatim}
m.BandSVD().Thresh(RT thresh)
m.BandSVD().Top(size_t nsing)
RT m.BandSVD().Condition()
size_t m.BandSVD().GetKMax()
\end{verbatim}

There are also the \tt{SVS}, \tt{SVU}, and \tt{SVV} options for the decomposition
which cannot be used for division, but which can access parts of the decomposition.

\end{enumerate}

The routines 
\begin{verbatim}
m.SaveDiv()
m.SetDiv()
m.ReSetDiv()
m.UnSetDiv()
m.DivideInPlace()
\end{verbatim}
work the same as for regular \tt{Matrix}es.

\subsubsection{Input/Output}

The simplest output is the usual:
\begin{verbatim}
os << m
\end{verbatim}
where \tt{os} is any \tt{std::ostream}.
The output format is the same as for a \tt{Matrix}.

There is also a compact format for a \tt{BandMatrix}:
\begin{verbatim}
m.WriteCompact(os)
\end{verbatim}
outputs in the format:
\begin{verbatim}
B nrows ncols nlo nhi
( m(0,0)  m(0,1)  m(0,2) ... m(0,nhi) )
( m(1,0)  m(1,1)  m(1,2) ... m(1,nhi+1) )
...
( m(nlo,0)  m(nlo,1) ...  m(nlo,nlo+nhi) )
...
( m(nrows-nhi-1,nrows-nlo-nhi-1) ... m(nrows-nhi-1,ncols-1) )
...
( m(nrows-1,nrows-nlo-1)  ... m(nrows-1,ncols-1) )
\end{verbatim}

The same compact format can be read back in the usual two ways:
\begin{verbatim}
tmv::BandMatrix<T> b(nrows,ncols,nlo,nhi);
is >> b;
std::auto_ptr<tmv::BandMatrix<T> > pb;
is >> pb;
\end{verbatim}

One can write small values as 0 with
\begin{verbatim}
m.Write(std::ostream& os, RT thresh)
m.WriteCompact(std::ostream& os, RT thresh)
\end{verbatim}

\subsection{Symmetric and hermitian band matrices}

The \tt{SymBandMatrix} class is our symmetric band matrix, which combines
the properties of SymMatrix and BandMatrix; it has a banded structure and
$m = m^T$.  Likewise \tt{HermBandMatrix} is our hermitian band matrix for
which $m = m^\dagger$.

Both classes inherit from \tt{GenSymBandMatrix}, which in turn inherits from 
\tt{BaseMatrix}.  The various views and composite classes described below 
also inherit from \tt{GenSymBandMatrix}.

As with the documentation for \tt{SymMatrix}/\tt{HermMatrix}, the descriptions
below will only be written for \tt{SymBandMatrix} with the implication
that a \tt{HermBandMatrix} has the same functionality, but with the calculations
appropriate for a hermitian matrix, rather than symmetric.

One general caveat about complex \tt{HermBandMatrix} calculations is that 
the diagonal elements should all be real.  Some calculations assume this, and 
only use the real part of the diagonal elements.  Other calculations use the 
full complex value as it is in memory.  Therefore, if you set a diagonal element
to a non-real value at some point, the results will likely be wrong in 
unpredictable ways.  Plus of course, your matrix will not actually be 
hermitian any more, so the right answer is undefined in any case.

All the \tt{SymBandMatrix} and \tt{HermBandMatrix} routines are included by:
\begin{verbatim}
#include "TMV_SymBand.h"
\end{verbatim}

In addition to the \tt{T} template parameter, there are three other template 
parameters: \tt{uplo} which can be either \tt{tmv::Upper} or \tt{tmv::Lower},
\tt{stor} which can be one of \tt{tmv::RowMajor}, \tt{tmv::ColMajor}, or 
\tt{tmv::DiagMajor}, and \tt{index} which can be either \tt{tmv::CStyle}
or \tt{tmv::FortranStyle}.  The default values for these are 
\tt{Upper}, \tt{RowMajor}, \tt{CStyle}.
The default values for these are \tt{RowMajor} and \tt{CStyle}.

The storage size required is the same as for the \tt{BandMatrix} of
the upper or lower band portion.
As with square band matrices, 
all three storage methods always need the same amount of memory, so the 
decision about which method to use should generally be based on performance 
considerations rather than memory usage.
The speed of the various matrix operations are different for the different storage 
types.  If the matrix calculation speed is important, it may be worth trying 
all three to see which is fastest for the operations you are using.

Also, as with \tt{BandMatrix}, the storage for \tt{Lower}, \tt{DiagMajor} 
does not start with the upper left element as usual.
Rather, it starts at the start of the lowest sub-diagonal.  
So for the constructors which 
take the matrix information from an array (\tt{T*} or \tt{vector<T>}), 
the start of the array needs to be at the start of the lowest sub-diagonal.  

\subsubsection{Constructors}

\begin{itemize}
\item 
\begin{verbatim}
tmv::SymBandMatrix<T,uplo,stor,index> m(size_t n, int nlo)
\end{verbatim}
Makes an \tt{n} $\times$ \tt{n} \tt{SymBandMatrix} with 
\tt{nlo} off-diagonals
and with \underline{uninitialized} values.
If debugging is turned on (i.e. not turned off 
with -DNDEBUG), then the values are initialized to 888.

\item
\begin{verbatim}
tmv::SymBandMatrix<T,uplo,stor,index> m(size_t n, int nlo, T x)
\end{verbatim}
Makes an \tt{n} $\times$ \tt{n} \tt{SymBandMatrix} \tt{nlo} off-diagonals
and with all values equal to \tt{x}.

\item
\begin{verbatim}
tmv::SymBandMatrix<T,uplo,stor,index> m(size_t n, int nlo, 
        const T* mm)
tmv::SymBandMatrix<T,uplo,stor,index> m(size_t n, int nlo, 
        const std::vector<T>& mm)
\end{verbatim}
Makes a \tt{SymBandMatrix} which copies the elements from \tt{mm}.  

\item 
\begin{verbatim}
tmv::SymBandMatrix<T,uplo,stor,index> m(const GenMatrix<T>& mm, 
        int nlo)
tmv::SymBandMatrix<T,uplo,stor,index> m(const GenSymMatrix<T>& mm, 
        int nlo)
tmv::SymBandMatrix<T,uplo,stor,index> m(const GenBandMatrix<T>& mm, 
        int nlo)
tmv::SymBandMatrix<T,uplo,stor,index> m(const GenSymBandMatrix<T>& mm, 
        int nlo)
\end{verbatim}
Makes an \tt{SymBandMatrix} which copies the corresponding values of \tt{mm}.  
For the last two, \tt{nlo} must not be larger than the number of upper
or lower bands in \tt{mm}.

\item
\begin{verbatim}
tmv::SymBandMatrixView<T> m = 
        SymBandMatrixViewOf(Matrix<T>& mm, UpLoType uplo, int nlo)
tmv::SymBandMatrixView<T> m = 
        SymBandMatrixViewOf(const MatrixView<T>& mm, UpLoType uplo, 
        int nlo)
tmv::ConstSymBandMatrixView<T> m = 
        SymBandMatrixViewOf(const GenMatrix<T>& m, UpLoType uplo, 
        int nlo)
tmv::SymBandMatrixView<T> m = 
        SymBandMatrixViewOf(SymMatrix<T>& mm, int nlo)
tmv::SymBandMatrixView<T> m = 
        SymBandMatrixViewOf(const SymMatrixView<T>& mm, int nlo)
tmv::ConstSymBandMatrixView<T> m = 
        SymBandMatrixViewOf(const GenSymMatrix<T>& m, int nlo)
tmv::SymBandMatrixView<T> m = 
        SymBandMatrixViewOf(BandMatrix<T>& mm, UpLoType uplo, int nlo)
tmv::SymBandMatrixView<T> m = 
        SymBandMatrixViewOf(const BandMatrixView<T>& mm, UpLoType uplo, 
        int nlo)
tmv::ConstSymBandMatrixView<T> m = 
        SymBandMatrixViewOf(const GenBandMatrix<T>& m, UpLoType uplo, 
        int nlo)
\end{verbatim}
Make an \tt{SymBandMatrixView} of the corresponding portion of \tt{mm}.  
To view these as a hermitian band matrix, use the command,
\tt{HermBandMatrixViewOf} instead.
For the views from \tt{BandMatrix}, the parameter \tt{nlo} may be 
omitted, in which case either \tt{m.nlo()} or \tt{m.nhi()} is used 
according to whether \tt{uplo} is \tt{Lower} or \tt{Upper} respectively.

\item
\begin{verbatim}
tmv::SymBandMatrixView<T> m = 
        tmv::SymBandMatrixViewOf(T* mm, size_t n, int nlo, 
          UpLoType uplo, StorageType stor)
tmv::ConstSymBandMatrixView<T> m = 
        tmv::SymBandMatrixViewOf(const T* mm, size_t n, int nlo, 
          UpLoType uplo, StorageType stor)
tmv::SymBandMatrixView<T> m = 
        tmv::HermBandMatrixViewOf(T* mm, size_t n, int nlo, 
          UpLoType uplo, StorageType stor)
tmv::ConstSymBandMatrixView<T> m = 
        tmv::HermBandMatrixViewOf(const T* mm, size_t n, int nlo, 
          UpLoType uplo, StorageType stor)
\end{verbatim}
Make a \tt{SymBandMatrixView} using the actual memory elements, \tt{mm}.

\item
\begin{verbatim}
tmv::SymBandMatrix<T,uplo,stor> m = SymTriDiagMatrix<uplo>(
        const GenVector<T>& v1, const GenVector<T>& v2)
tmv::SymBandMatrix<T,uplo,stor> m = SymTriDiagMatrix<uplo>(
        const GenVector<T>& v1, const GenVector<T>& v2)
\end{verbatim}
Shorthand to create a symmetric tri-diagonal \tt{BandMatrix} 
if you already have the \tt{Vector}s.  
The main diagonal is \tt{v1} and the off-diagonal is \tt{v2}.

\item
\begin{verbatim}
tmv::SymBandMatrix<T> m1 = m2
tmv::SymBandMatrix<T> m1(const GenSymBandMatrix<T2>& m2)
\end{verbatim}
Copy the \tt{SymBandMatrix m2}, which may be of any type \tt{T2} so long
as values of type \tt{T2} are convertible into type \tt{T}.
The assignment operator has the same flexibility.

\end{itemize}

\subsubsection{Access}

\begin{verbatim}
m.nrows() = m.ncols() = m.colsize() = m.rowsize() = m.size()
m.nlo() = m.nhi()
m(i,j)
m.row(size_t i, size_t j1, size_t j2)
m.col(size_t i, size_t j1, size_t j2)
m.diag()
m.diag(int i)
m.diag(int i, size_t k1, size_t k2)
\end{verbatim}
Again, the versions of \tt{row} and \tt{col} with only one argument are
missing, since the full row or column isn't accessible as a \tt{VectorView}.
You must specify a valid range within the row or column that you want, 
given the banded storage of \tt{m}.

\begin{verbatim}
m.SubVector(size_t i, size_t j, int istep, int jstep, size_t size)
m.SubMatrix(int i1, int i2, int j1, int j2)
m.SubMatrix(int i1, int i2, int j1, int j2, int istep, int jstep)
m.SubBandMatrix(int i1, int i2, int j1, int j2, int newnlo, int newhi)
m.SubBandMatrix(int i1, int i2, int j1, int j2, int newnlo, int newhi, 
        int istep, int jstep)
m.Diags(int k1, int k2)
m.UpperBand()
m.LowerBand()
\end{verbatim}
These work the same as for a \tt{BandMatrix}, except that the entire
subvector or submatrix must be completely within the upper or lower band.

\begin{verbatim}
m.SubSymMatrix(int i1, int i2)
m.SubSymMatrix(int i1, int i2, int istep)
m.SubSymBandMatrix(int i1, int i2, int newnlo=m.nlo())
m.SubSymBandMatrix(int i1, int i2, int newnlo, int istep)
\end{verbatim}
These return a view of the \tt{SymMatrix} or \tt{SymBandMatrix} which runs
from \tt{i1} to \tt{i2} along the diagonal with an optional step,
and includes the off-diagonals in the same rows/cols.  For the first two,
the \tt{SymMatrix} must be completely with the band.

\begin{verbatim}
m.UpperBandOff()
m.LowerBandOff()
\end{verbatim}
These return a \tt{BandMatrixView} of only the off-diagonals of either the
upper or lower half of the matrix.  With a \tt{SymMatrix} some algorithms
use the combination \tt{m.UpperTri().OffDiag()}.  The corresponding
operation for a \tt{SymBandMatrix} would be \tt{m.UpperBand().OffDiag()},
but a \tt{BandMatrix} does not have the method \tt{OffDiag}.  
Hence, these methods were added to do this operation directly.

\begin{verbatim}
m.Real()
m.Imag()
m.Transpose()
m.Conjugate()
m.Adjoint()
m.View()
\end{verbatim}
Note that the imaginary part of a complex hermitian band matrix is
anti-symmetric, so \tt{m.Imag()} is illegal for a \tt{HermBandMatrix}.
If you need to manipulate the imaginary part of a \tt{HermMatrix}, 
you could use
\tt{m.UpperBandOff().Imag()} 
(since all the diagonal elements are real).

\subsubsection{Functions}

\begin{verbatim}
RT m.Norm1() = Norm1(m)
RT m.Norm2() = Norm2(m)
RT m.NormInf() = NormInf(m)
RT m.NormF() = NormF(m) = m.Norm() = Norm(m)
RT m.NormSq() = NormSq(m)
RT m.MaxAbsElement() = MaxAbsElement(m)
T m.Trace() = Trace(m)
T m.Det() = Det(m)
minv = m.Inverse() = Inverse(m)
m.Inverse(Matrix<T>& minv)
m.Inverse(SymMatrix<T>& minv)
m.InverseATA(Matrix<T>& cov)
\end{verbatim}
The inverse of a \tt{SymBandMatrix} is not (in general) banded.  
However, it is symmetric (or hermitian).
So \tt{minv} here must be a regular \tt{Matrix} or \tt{SymMatrix}.

\begin{verbatim}
m.Zero()
m.SetAllTo(T x)
m.Clip(RT thresh)
m.SetToIdentity(T x = 1)
m.ConjugateSelf()
m.TransposeSelf()
Swap(m1,m2)
\end{verbatim}

\subsubsection{Arithmetic}

In addition to \tt{x}, \tt{v}, \tt{m}, \tt{b} and \tt{s} from before, 
we now add \tt{sb} for a \tt{SymBandMatrix}.

\begin{verbatim}
sb2 = -sb1

sb2 = x * sb1
sb2 = sb1 [*/] x

sb3 = sb1 [+-] sb2
m2 = m1 [+-] sb
m2 = sb [+-] m1
b2 = b1 [+-] sb
b2 = sb [+-] b1
s2 = s1 [+-] sb
s2 = sb [+-] s1

sb [*/]= x

sb2 [+-]= sb1
m [+-]= sb
b [+-]= sb
s [+-]= sb

v2 = sb * v1
v2 = v1 * sb
v *= sb

b = sb1 * sb2
m2 = sb * m1
m2 = m1 * sb
m *= sb
b2 = sb * b1
b2 = b1 * sb
b *= sb
m2 = sb * s1
m2 = s1 * sb

sb2 = sb1 [+-] x
sb2 = x [+-] sb1
sb [+-]= x
\end{verbatim}

\subsubsection{Division}

The division operations are:
\begin{verbatim}
v2 = v1 [/%] b
m2 = m1 [/%] b
m2 = b [/%] m1
m = b1 [/%] b2
v [/%]= b
m [/%]= b
\end{verbatim}

\tt{SymBandMatrix} has three possible choices for the division decomposition:
\begin{enumerate}
\item
\tt{m.DivideUsing(tmv::LU)} actually does the \tt{BandMatrix} version of 
LU, rather than a Bunch-Kauffman algorithm like for \tt{SymMatrix}.  The
reason is that the pivots in the Bunch-Kauffman algorithm can arbitrarily
expand the band width required to hold the information.  The generic
banded LU algorithm is limited to 3*\tt{nlo}+1 bands.

To access this decomposition, use:
\begin{verbatim}
bool m.LUD().IsTrans()
tmv::LowerTriMatrix<T,UnitDiag> m.LUD().GetL()
tmv::ConstBandMatrixView<T> m.LUD().GetU()
size_t* m.LUD().GetP()
\end{verbatim}
The following should result in a matrix numerically very close to \tt{m}.
\begin{verbatim}
tmv::Matrix<T> m2(m.nrows,m.ncols);
tmv::MatrixView<T> m2v = 
        m.LUD().IsTrans() ? m2.Transpose() : m2.View();
m2v = m.LUD().GetL() * m.LUD().GetU();
m2v.ReversePermuteRows(m.LUD().GetP());
\end{verbatim}

\item
\tt{m.DivideUsing(tmv::CH)} will perform a Cholesky decomposition.  
\tt{m} must be hermitian (or real symmetric) to use \tt{CH}, since that is the
only kind of matrix that has a Cholesky decomposition.  

As with a regular \tt{SymMatrix},
the only real advantage of Cholesky over Bunch-Kauffman is speed.  If you know your 
matrix is positive-definite, the Cholesky decomposition is the fastest way to 
do division.

If \tt{m} is tri-diagonal (i.e. \tt{nlo} = 1), then we use a slightly 
different algorithm, which avoids the square roots usually required for 
Cholesky.  Namely, we form the decomposition $m = LDL^\dagger$, where $L$ is a
unit-diagonal lower banded matrix with 1 sub-diagonal, and $D$ is diagonal.

If \tt{m} has \tt{nlo} $> 1$, then we just use a normal Cholesky algorithm
where $m = LL^\dagger$ and $L$ is lower banded with the same \tt{nlo} as
\tt{m}.

Both versions of the algorithm are accessed with the same methods:
\begin{verbatim}
BandMatrix<T> m.CHD().GetL()
DiagMatrix<T> m.CHD().GetD()
\end{verbatim}
with $L$ being made unit-diagonal or $D$ being set to the identity matrix
as appropriate.

The following should result in a matrix numerically very close to \tt{m}.
\begin{verbatim}
Matrix<T> m2 = m.CHD().GetL() * m.CHD().GetD() * 
        m.CHD().GetL().Adjoint()
\end{verbatim}

\item
\tt{m.DivideUsing(tmv::SV)} will perform a singular value decomposition.

As with regular \tt{SymMatrix}es, the singular value decomposition of 
hermitian band matrices is also an eigenvalue decomposition: $m = USU^\dagger$.
Thus, $S$ may have negative values, and is therefore not quite a real
SVD (for which all the singular values are require to be positive).

To access this decomposition, use:
\begin{verbatim}
ConstMatrixView<T> m.SVD().GetU()
ConstVectorView<RT> m.SVD().GetS()
ConstMatrixView<T> m.SVD().GetV()
\end{verbatim}
The following should result in a matrix numerically very close to \tt{m}.
\begin{verbatim}
Matrix<T> m2 = m.SVD().GetU() * DiagMatrixViewOf(m.SVD().GetS()) * 
        m.SVD().GetV()
\end{verbatim}

For a complex, symmetric, banded \tt{m}, the accessor is through
\tt{m.SymSVD()}, just like for a complex \tt{SymMatrix} object.

Both versions also have the same control and access routines as a regular SVD:
\begin{verbatim}
m.SVD().Thresh(RT thresh)
m.SVD().Top(size_t nsing)
RT m.SVD().Condition()
size_t m.SVD().GetKMax()
\end{verbatim}
(Likewise for m.SymSVD().)

There are also the \tt{SVS}, \tt{SVU}, and \tt{SVV} options for the decomposition
which cannot be used for division, but which can access parts of the decomposition.

\end{enumerate}

For all three \tt{DivType}s, the routines 
\begin{verbatim}
m.SaveDiv()
m.SetDiv()
m.ReSetDiv()
m.UnSetDiv()
m.DivideInPlace()
\end{verbatim}
work the same as for regular \tt{Matrix}es.

\subsubsection{Input/Output}

The simplest output is the usual:
\begin{verbatim}
os << m
\end{verbatim}
where \tt{os} is any \tt{std::ostream}.
The output format is the same as for a \tt{Matrix}.

There is also a compact format for a \tt{BandMatrix}:
\begin{verbatim}
m.WriteCompact(os)
\end{verbatim}
outputs in the format:
\begin{verbatim}
hB/sB n nlo
( m(0,0)  )
( m(1,0)  m(1,1)  )
...
( m(nlo,0)  m(nlo,1) ...  m(nlo,nlo) )
( m(nlo+1,1)  m(nlo+1,2) ...  m(nlo+1,nlo+1) )
...
( m(n-1,n-nlo-1)  ... m(n-1,n-1) )
\end{verbatim}
where \tt{hB/sB} means \underline{either} \tt{hB} or \tt{sB}, 
which indicates whether the matrix is hermitian or symmetric.

The same compact format can be read back in the usual two ways:
\begin{verbatim}
tmv::SymBandMatrix<T> s(n,nlo);
tmv::HermBandMatrix<T> h(n,nlo);
is >> s >> h;
std::auto_ptr<tmv::SymBandMatrix<T> > ps;
std::auto_ptr<tmv::HermBandMatrix<T> > ph;
is >> ps >> ph;
\end{verbatim}

One can write small values as 0 with
\begin{verbatim}
m.Write(std::ostream& os, RT thresh)
m.WriteCompact(std::ostream& os, RT thresh)
\end{verbatim}

\newpage
\section{Errors and Exceptions}

There are two kinds of errors that I looks for.  The first are 
coding errors.  Some examples are:
\begin{itemize}
\item
Trying to access elements outside the range of a \tt{Vector} or \tt{Matrix}.
\item
Trying to add to \tt{Vector}s or \tt{Matrix}es which are different sizes.
\item
Trying to multiply a \tt{Matrix} by a \tt{Vector} where the number of columns 
in the \tt{Matrix} doesn't match the size of the \tt{Vector}.
\item
Viewing a \tt{Matrix} as a \tt{HermMatrix} when the diagonal isn't real.
\item
Calling \tt{m.LUD()}, \tt{m.QRD()}, etc. for a \tt{Matrix} which does not have that 
decomposition set (and saved).
\end{itemize}

Basically, these are all things which ideally the compiler would catch, since they 
are programming errors, but either they are impossible/impractical to implement
in this way, or I just haven't bothered to do so.

I check for all of these (and similar) errors using assert statements.
If these asserts fail, it should mean that
the programmer made a mistake in the code.  (Unless I've made a 
mistake in the TMV code, of course.)

Once the code is working, you can make the code slightly faster by 
compiling with -DNDEBUG.  I say slightly, since most of these checks 
are pretty innocuous.  And most of the computing time is usually in the depths
of the various algorithms, not in these $O(1)$ time checks that the dimensions
are correct and such.  

The other kind of error checked for by the code is 
where the data don't behave in the way
the programmer expected.  Here is a (complete) list of these errors:
\begin{itemize}
\item
A singular matrix is encountered in a division routine that cannot handle it.
\item
An input file has the wrong format.
\item
A Cholesky decomposition is attempted for a
hermitian matrix which isn't positive definite.
\end{itemize}

These errors are always checked for even if -DNDEBUG is used.
That's because they are not problems in the code per se, but rather
problems with the data or files used by the code.  So they could
still happen even after the code has been thoroughly tested.

All errors in the TMV library are indicated by throwing an object of type
\tt{tmv::Error}.  If you decide to catch it, you can determine what went
wrong by printing it:
\begin{verbatim}
catch (tmv::Error& e) {
  std::cerr << e << std::endl;
}
\end{verbatim}
If you catch the error by value (i.e. without the \tt{\&}), it will print out
a single line description.  If you catch it by reference (as above), it may 
print out more information about the problem.

Also, \tt{tmv::Error} derives from \tt{std::exception} and overrides the \tt{what()}
method, so any program which catches these will catch \tt{tmv::Error} as well.

If you want to be more specific, there are a number of classes which derive from
\tt{Error}:
\begin{itemize}
\item
\tt{tmv::FailedAssert} indicates that one of the assert statements failed.
Since these are coding errors,
if you catch this one, you'll probably just want to print out the error and abort
the program so you can fix the bug.  In addition to printing the text of the
assert statement that failed (if you catch by reference), it will also indicate
the file and line number as normal assert macros do.  Unfortunately, it gives the 
line number in the TMV code, rather than in your own code, but hopefully seeing
which function in TMV found the problem will help you figure out which line
in your own code was incorrect.

If you are absolutely sure that the assert failed due to a bug in the TMV code
rather than your own code, please let me know at \myemail.

\item
\tt{tmv::Singular} indicates that you tried to invert or divide by a matrix
which is (numerically) singular.  This may be useful to catch specifically,
since you may want to do something different when you encounter a singular
matrix.  Note however, that this only detects \underline{exactly} singular
matrices.  If a matrix is numerically close to singular, but no actual 
zeros are found, then no error will be thrown, but your results will be 
unreliable.

\item
\tt{tmv::ReadError} indicates that there was some problem reading in a matrix or 
vector from an \tt{istream} input.  If you catch this by reference and write it,
it will give you a fairly specific description of what the problem was as well as
writing the part of the matrix or vector which was read in successfully.

\item
\tt{tmv::NonPosDef} indicates that you tried to do a Cholesky decomposition of
a hermitian matrix which was not positive definite.  I suspect that this is the 
most useful exception to catch specifically, as opposed to just via \tt{tmv::Error}.

For example, the fastest algorithm for determining whether a matrix is 
(at least numerically) positive definite is to try the Cholesky decomposition
and catch this exception.  To wit:
\begin{verbatim}
bool IsPosDef(const tmv::GenSymMatrix<T>& m)
{
  assert(m.isherm());
  try {
    tmv::ConstSymMatrixView<T> mview = m;
    mview.DivideUsing(tmv::CH);
    mview.SetDiv();
  }
  catch (tmv::NonPosDef) {
    return false;
  }
  return true;
}
\end{verbatim}
(Using the \tt{ConstSymMatrixView} is a trick to avoid destroying any decomposition
that is already set for \tt{m}.)

Or you might want to use Cholesky when possible and Bunch-Kauffman otherwise:
\begin{verbatim}
try {
  m.DivideUsing(tmv::CH);
  m.SetDiv();
}
catch (tmv::NonPosDef) {
  m.DivideUsing(tmv::LU);
  m.SetDiv();
}
x = b/m;
\end{verbatim}

Note, however, that the speed difference between the two algorithms is only about
20\% for typical matrices.  So if a significant fraction of your matrices are 
not positive definite, you are probably better off always using LU.  

Code like that given above would probably only be useful when all of your 
matrices \underline{should} be positive definite in exact arithmetic, but you 
want to guard against one failing the Cholesky decomposition due to round-off
errors.

It is also worth mentioning that the routine \tt{QR\_Downdate} described in \S\ref{downdate}
below will also throw the exception \tt{NonPosDef} when it fails.

\end{itemize}

\newpage
\section{Advanced Usage}

\subsection{Element-by-element product}

The two usual kinds of multiplication for vectors are the inner product and 
the outer product, which result in a scalar and a matrix respectively.
However, occasionally you may want to multiply each element in a vector by the
corresponding element in another vector: $v(i) = v(i) * w(i)$.

There are two functions that should provide all of this kind of functionality
for you:
\begin{verbatim}
ElementProd(T x, const GenVector<T1>& v1, const VectorView<T>& v2);
AddElementProd(T x, const GenVector<T1>& v1, const GenVector<T2>& v2,
        const VectorView<T>& v3)
\end{verbatim}
The first performs $v2(i) = x * v1(i) * v2(i)$, and the second performs
$v3(i) = v3(i) + x * v1(i) * v2(i)$ for $i = 0 ... (N-1)$ (where $N$ is the 
size of the vectors).

There is no operator overloading for \tt{Vector}s that would be equivalent to 
these expressions.
But they are actually equivalent to the following forms:
\begin{verbatim}
v2 *= x * DiagMatrixViewOf(v1);
v3 += x * DiagMatrixViewOf(v1) * v2;
\end{verbatim}
respectively.  In fact, these statements inline to the above function calls
automatically.  Depending on you preference and the meanings of your vectors,
these statements may or may not be clearer as to what you are doing.

There are also corresponding functions for \tt{Matrix} and for each of the special
matrix types:
\begin{verbatim}
ElementProd(T x, const GenMatrix<T1>& m1, const MatrixView<T>& m2);
AddElementProd(T x, const GenMatrix<T1>& m1, const GenMatrix<T2>& m2,
        const MatrixView<T>& m3);
\end{verbatim}
Likewise for the other special matrix classes.  The first performs 
$m2(i,j) = x * m1(i,j) * m2(i,j)$, and the second performs
$m3(i,j) = m3(i,j) + x * m1(i,j) * m2(i,j)$ for every $i,j$ in the matrix.

These don't have any \tt{DiagMatrixViewOf} version, since the corresponding 
concept would require a four-dimensional tensor, and (so far!) the TMV library
just deals with one- and two-dimensional objects.

The matrices all have to be the same size and shape, but can have any 
(i.e. not necessarily the same) storage method.  However, the routines are faster
if the matrices all use the same storage.

\subsection{QR decomposition and update}

Often, it can be useful to create and deal with the QR decomposition directly,
rather than 
just relying on the division routines.  Primarily, this is because of the possibility
of updating or ``downdating'' the resulting $R$ matrix.  

If you are doing a 
least-square fit to a large number of linear equations, the matrix version of this
is: $A x = b$, where $A$ is a matrix with more rows than columns (i.e. it is tall).
($A$ is commonly called the ``design matrix''.)
It may be the case that you
have more rows (i.e. constraints) than would allow the entire matrix to fit in memory.  

In this case it may be tempting to use the so-called normal equation instead: 
\begin{align*}
A^\dagger A x &= A^\dagger b \\
x & = (A^\dagger A)^{-1} A^\dagger b
\end{align*}
This equation theoretically gives the same 
solution as using the QR decomposition on the original design matrix.
However, it can be shown that the condition of $A^\dagger A$ is the 
\underline{square} of the condition of $A$.  Since larger condition values
lead to larger numerical instabilities and round-off problems, a mildly
ill-conditioned matrix is made much worse by this procedure.

A better solution is to use the QR decomposition, $A = QR$, so
\begin{align*}
Q R x &= b \\
x &= R^{-1} Q^\dagger b
\end{align*}
But if $A$ is too large to fit in memory, then so is $Q$.

A compromise solution, which is not quite as good as doing the full QR decomposition,
but is better than using the normal equation, is to just calculate the $R$ of the
QR decomposition, and not $Q$.  Then:
\begin{align*}
A^\dagger A x &= A^\dagger b \\
R^\dagger Q^\dagger Q R x = R^\dagger R x &= A^\dagger b \\
x &= R^{-1} (R^\dagger)^{-1} A^\dagger b
\end{align*}

Calculating $R$ directly from $A$ is numerically much more stable than 
calculating it through a Cholesky decomposition of $A^\dagger A$.
So this method produces a more accurate answer for $x$ than the normal equation does.

How do you calculate $R$?

First, we use the lemma that a product of unitary matrices is also unitary.  
This implies that if we can calculate
something like: $A = Q_0 Q_1 Q_2 ... Q_n R$, then this $R$ is the one we want.

So, consider breaking $A$ into a submatrix, $A_0$, which can fit into memory, 
plus the remainder, $A_1$.
\begin{equation*}
A = \left(\begin{array}{c}A_0 \\A_1\end{array}\right)
\end{equation*}
First perform a QR decomposition of $A_0 = Q_0 R_0$.  Then we have:
\begin{align*}
A &= \left(\begin{array}{c}Q_0 R_0 \\ A_1 \end{array}\right) \\
&= \left(\begin{array}{cc}Q_0 & 0 \\ 0 & 1\end{array}\right) 
      \left(\begin{array}{c}R_0 \\ A_1 \end{array}\right) \\
&\equiv Q_0^\prime A_1^\prime
\end{align*}

Assuming that $A_0$ has more rows than columns, 
then $A_1^\prime$ has fewer rows than
the original matrix $A$.  So we can iterate this process until the 
resulting matrix can fit in memory, and we can perform the final QR update
to get the final value of $R$.

For numerical reasons, 
the fewer such iterations you do, the better.  So you should try to make the 
top submatrix as large as possible each time, given the amount of memory
available.

The commands to do this with the TMV library are:
\begin{verbatim}
QR_Decompose(const MatrixView<T>& A0, const VectorView<T>& Qbeta);
UpperTriMatrix<T> R = UpperTriMatrixViewOf(A0);
QR_Update(const UpperTriMatrixView<T>& R, const GenMatrix<T>& A1)
QR_Update(const UpperTriMatrixView<T>& R, const GenMatrix<T>& A2)
...
\end{verbatim}

The first subroutine performs the QR decomposition of \tt{A} in place, with $Q$ 
stored in a compact form below the diagonal of \tt{A} and in the extra
vector \tt{Qbeta}.  Since we don't need $Q$, you can throw away
\tt{Qbeta} after this call.  \tt{R} is returned in the upper triangle of \tt{A}
as shown.

Then subsequent updates would use the second subroutine, \tt{QR\_Update(R,An)}.
For this routine, \tt{A1}, \tt{A2}, etc. are the additional rows to 
add to the QR decomposition in each step.  

The solution equation, written above, also needs the quantity $A^\dagger b$, which
can be accumulated in the same blocks:
\begin{equation}
A^\dagger b = A_1^\dagger b_1 + A_2^\dagger b_2 + ....
\end{equation}
This, combined with the calculation of $R$, allows us to determine $x$ using the above formula.

\subsection{QR downdate}
\label{downdate}

When performing a least-square solution to some model as described in the 
previous section, 
it is common to do some kind of outlier rejection to remove data which
seem not to be applicable to the model - things like spurious measurements
and such.
For this, we basically want the opposite of a QR update - instead we want to 
find the QR decomposition that results from
removing a few rows from $A$.  This is called a QR ``downdate'', and is performed
using the subroutine:
\begin{verbatim}
void QR_Downdate(const UpperTriMatrixView<T>& R, const GenMatrix<T>& M)
\end{verbatim}
Now \tt{M} represents the rows from the original matrix to remove from the 
QR decomposition.

It is possible for the downdate to fail if the matrix $M$ does not represent rows
of the matrix which was originally used to create $R$.
Furthermore,
with round-off errors, the error may still result with actual rows from the 
original $A$
if $R$ gets too close to singular.  In this case, \tt{QR\_Downdate(R,M)} throws
a \tt{NonPosDef} exception.  This might seem like a strange choice, but the 
logic is that $R^\dagger R$ is the Cholesky decomposition of $A^\dagger A$,
and \tt{QR\_Downdate(R,M)} basically updates $R$ to be the Cholesky decomposition
of $A^\dagger A - M^\dagger M$.  The procedure fails (and throws) when this latter 
matrix is found not to be positive definite.

\subsection{Other SymMatrix operations}

There are three arithmetic routines which we provide for \tt{SymMatrix},
which do not have
any corresponding shorthand with the usual arithmetic operators.

The first two are:
\begin{verbatim}
tmv::Rank2Update<bool add>(T x1, const GenVector<T1>& v1, 
        const GenVector<T2>& v2, const SymMatrixView<T>& m)
tmv::Rank2KUpdate<bool add>(T x1, const GenMatrix<T1>& m1, 
        const GenMatrix<T2>& m2, const SymMatrixView<T>& m3)
\end{verbatim}
They are similar to the \tt{Rank1Update} and \tt{RankKUpdate} routines
which are implemented via the \texttt{\^} operator.

A rank-2 update calculates
\begin{verbatim}
m (+=) x1 * ((v1 ^ v2) + (v2 ^ v1))} 
m (+=) x1 * (v1 ^ v2.Conjugate()) + conj(x1) * (v2 ^ v1.Conjugate())}
\end{verbatim}
for a symmetric or hermitian \tt{m} respectively,
where ``(+=)'' means ``+='' if \tt{add} is \tt{true} and ``='' 
if \tt{add} is \tt{false}.
Likewise, a rank-2k update calculates:
\begin{verbatim}
m3 (+=) x1 * (m1 * m2.Transpose() + m2 * m1.Transpose())
m3 (+=) x1 * m1 * m2.Adjoint() + conj(x1) * m2 * m1.Adjoint()
\end{verbatim}
for a symmetric or hermitian \tt{m} respectively.

We don't have an arithmetic operator 
shorthand for these, because, as you can see, the operator
overloading required would be quite complicated.  
And since they are pretty rare, I decided to just let the programmer 
call the routines explicitly.

The other routine is:
\begin{verbatim}
tmv::SymMultMM<bool add>(T x1, const GenMatrix<T>& m1, 
        const GenMatrix<T>& m2, const SymMatrixView<T>& m3)
\end{verbatim}
This calculates the usual generalized matrix product:
\tt{m3 (+=) x1*m1*m2}, but it basically
asserts that the product \tt{m1*m2} is symmetric (or hermitian as appropriate).

Since a matrix product is not in general symmetric, I decided not to allow 
this operation with just the usual operators to prevent the user from doing 
this accidentally.  However, there are times when the 
programmer can know that the product should be (at least numerically close to)
symmetric and that this calculation is ok.  Therefore it is provided as a subroutine.

\subsection{Iterators}

We mentioned that the iterators through a \tt{Vector} are:
\begin{verbatim}
typename tmv::Vector<T>::iterator
typename tmv::Vector<T>::const_iterator
typename tmv::Vector<T>::reverse_iterator
typename tmv::Vector<T>::const_reverse_iterator
\end{verbatim}
just like for standard library containers.  The specific types to which these
typedefs refer are:
\begin{verbatim}
tmv::VIt<T,tmv::Unit,tmv::NonConj>
tmv::CVIt<T,tmv::Unit,tmv::NonConj>
tmv::VIt<T,tmv::Step,tmv::NonConj>
tmv::CVIt<T,tmv::Step,tmv::NonConj>
\end{verbatim}
respectively.

\tt{VIt} is a mutable-iterator, and \tt{CVIt} is a const-iterator.  \tt{Unit} 
indicates that the step size is 1, while \tt{Step} allows for any step size
between successive elements (and is therefore slower).  For the reverse
iterators, the step size is -1.

This can be worth knowing if you are going to be optimizing code that uses
iterators of \tt{VectorView}s.
This is because their iterators are instead:
\begin{verbatim}
tmv::VIter<T>
tmv::CVIter<T>
\end{verbatim}
which always check the step size (rather than assuming unit steps) and always
keep track of a possible conjugation.

If you know that you are dealing with a view which is not conjugated, you can 
convert your iterator into one of the above \tt{VIt} or \tt{CVIt} types which will be 
faster, since they won't check the conjugation bit each time. 

Likewise, if you
know that it {\em is} conjugated, then you can use \tt{tmv::Conj} for the 
third template parameter above.  This indicates that the vector view really
refers to the conjugates of the values stored in the actual memory locations.

Also, if you know that your view has unit steps between elements, converting to 
an iterator with \tt{tmv::Unit} will iterate faster.  It is often faster to check
the step size once at the beginning of the routine and convert to a unit-step
iterator if possible.

All of these conversions can be done with a simple cast or constructor, such as:
\begin{verbatim}
if (v.step() == 1)
  for(VIt<float,Unit,NonConj> it = v.begin(); it != v.end(); it++)
    (*it) = sqrt(*it);
else
  for(VIt<float,Step,NonConj> it = v.begin(); it != v.end(); it++)
    (*it) = sqrt(*it);
\end{verbatim}

Regular \tt{Vector}s are always \tt{Unit} and \tt{NonConj}, so those iterators
are already fast without using the specific \tt{VIt} names. 
That is, you can just use \tt{Vector<T>::iterator} rather than \tt{VIt<T,Unit,NonConj>}
without any performance drop.

\subsection{Direct memory access}

We provide methods for accessing the memory of a matrix or vector directly.
This is especially useful for meshing the TMV objects with other libraries
(such as BLAS or LAPACK).  But it can also be useful for writing some
optimized code for a particular function.  

The pointer to the start of the memory for a vector can be obtained by:
\begin{verbatim}
T* v.ptr()
const T* v.cptr() const
\end{verbatim}

Using the direct memory access
requires that you know the spacing of the elements in memory and
(for views) whether the view is conjugated or not.  So we also provide:
\begin{verbatim}
int v.step() const
bool v.isconj() const
\end{verbatim}

For matrices, the corresponding routines actually return the upper-left element
of the matrix.  For some matrices, (e.g. \tt{DiagMajor BandMatrix}es) 
this is not actually the first element in memory.  We also need to know the 
step size in both directions:
\begin{verbatim}
T* m.ptr()
const T* m.cptr() const
int m.stepi() const
int m.stepj() const
bool m.isconj() const
bool m.isrm() const
bool m.iscm() const
\end{verbatim}
The step in the ``down'' direction along a column is \tt{stepi}, and the step to 
the ``right'' along a row is \tt{stepj}.
The last two check if a matrix is \tt{RowMajor} or \tt{ColMajor} respectively.

For band matrices, there are also:
\begin{verbatim}
int m.diagstep() const
bool m.isdm() const
\end{verbatim}
which return the step along the diagonal and whether the matrix is \tt{DiagMajor}.

For symmetric/hermitian matrices, there are some more methods:
\begin{verbatim}
bool m.isherm()
bool m.issym()
bool m.isupper()
\end{verbatim}
The first two both return \tt{true} for real symmetric matrices, but 
differentiate between symmetric and hermitian varieties for complex types.
The last one tells you whether the actual elements to be accessed are stored
in the upper triangle half of the matrix (true) or the lower (false).

\subsection{\tt{BaseMatrix} views}

If you are dealing with objects which are only known to be \tt{BaseMatrix}es
(i.e. they could be a \tt{Matrix} or a \tt{DiagMatrix} or a \tt{SymMatrix}, etc.),
then methods like \tt{m.Transpose()}, \tt{m.View()}, and such
can't know what kind of object to return.
So these methods can't be defined for a \tt{BaseMatrix}.  

Instead, we have the following virtual methods, 
which are available to a \tt{BaseMatrix}
object and are defined in each specific kind of matrix to return a pointer
to the right kind of object:
\begin{verbatim}
std::auto_ptr<tmv::BaseMatrix<T> > m.NewCopy()
std::auto_ptr<tmv::BaseMatrix<T> > m.NewView()
std::auto_ptr<tmv::BaseMatrix<T> > m.NewTranspose()
std::auto_ptr<tmv::BaseMatrix<T> > m.NewConjugate()
std::auto_ptr<tmv::BaseMatrix<T> > m.NewAdjoint()
std::auto_ptr<tmv::BaseMatrix<T> > m.NewInverse()
\end{verbatim}
\tt{NewCopy} and \tt{NewInverse} create new storage to store a copy of the 
matrix or its inverse, respectively.  The other four just return views of the current 
matrix.

\subsection{``Linear'' views}

Our matrices generally store the data contiguously in memory with all of the 
methods like \tt{row} and \tt{col} returning the appropriate slice through the
data.  Occasionally, though, it can be useful to treat the whole matrix
as a single vector of elements.  We use this internally for implementing routines
like \tt{SetAllTo} and matrix addition, among others.  These are faster than
accessing the data in ways that use the actual matrix structure.

Anyway, this same kind of access may be useful for some users of the library, 
so the following two methods are available:
\begin{verbatim}
tmv::VectorView<T> m.LinearView()
tmv::ConstVectorView<T> m.ConstLinearView()
\end{verbatim}
These return a view to the elements of a \tt{Matrix} as a single vector.  
It is always allowed for an actual \tt{Matrix}.  For a \tt{MatrixView} 
(or \tt{ConstMatrixView}), it is only allowed if all of the elements in the 
view are in one contiguous block of memory.  The helper function 
\tt{m.CanLinearize()} returns whether or not the above methods are legal.

The same methods are also defined for \tt{BandMatrix} (and corresponding views).
In this case, there are a few elements in memory which are not necessarily
defined, since they lie outside of the actual band structure, so some care
should be used depending on the application of the returned vector views.

The triangular and symmetric matrices have too much memory which is not
actually used by the matrix for these to be very useful, so we do not provide them.
When we eventually implement the packed storage varieties, these methods will
be provided for them.

A similar capability for \tt{Vector}s is:
\begin{verbatim}
tmv::VectorView<RT> v.Flatten()
\end{verbatim}
which returns a real view to the real and imaginary elements of a complex \tt{Vector},
which is required to have unit step.  The returned view has twice the 
length of \tt{v} and also has unit step.

This probably isn't very useful for most users either, but it is useful internally,
since it allows code such as:
\begin{verbatim}
tmv::Vector<complex<double> > v(10,complex<double>(1,3));
v *= 2.3;
\end{verbatim}
to call the BLAS routine \tt{dscal} with \tt{x=2.3}, rather than \tt{zscal}
with \tt{x=complex}<\tt{double}>\tt{(2.3,0.0)}, which would be slower.

\newpage
\section{Obtaining and Compiling the Library}

This code is licensed using the Gnu General Public License.  See 
\S\ref{license} below for more details.

The following are step-by-step instructions on how to obtain the code for the library,
set it up for your system, and compile it:
\begin{enumerate}

\item
Go to \website\ for a link to a tarball with all of the source code, and copy
it to the directory where you want to put the TMV library.

\item
Unpack the tarball:\\
\tt{gunzip tmv\tmvversion .tar.gz}\\
\tt{tar xf tmv\tmvversion .tar}

This will make a directory called tmv\tmvversion\ with the subdirectories:
\tt{doc}, \tt{examples}, \tt{include}, \tt{lib}, \tt{src}, \tt{test} 
along with the files \tt{INSTALL}, \tt{README} and \tt{Makefile} in the top 
directory.

\item
Edit the makefile.

The start of the makefile lists 4 things to specify: the compiler,
the include directories, the other flags to send to the compiler, 
and any necessary BLAS/LAPACK linkage flags.
The default setup is:
\begin{verbatim}
CC= g++
INCLUDE= -Iinclude
CFLAGS= $(INCLUDE) -O -DNOBLAS -DNDEBUG
BLASLIBS=
\end{verbatim}
but you will probably want to change this.

This default setup will compile using g++ without any BLAS or LAPACK library and 
with (a minimum level of) debugging.  It will place the library, \tt{libtmv.a},
and all the object files in the directory \tt{gcc-debug}.  And the test suite
executable will be placed in the current directory (where the makefile and code are).
This setup should work on any system with gcc, although it almost certainly 
won't be as fast as using an optimized BLAS library and/or a LAPACK library.

After these lines, there are several sets for different systems using various BLAS
and LAPACK versions, showcasing some of the compiler options described below
and giving examples of what you need for several common (or at least representative) systems.
If you have a system similar to one of these, then it should be a good 
starting point for you to figure out what you want to use.

Here are some compiler flags to consider using:
\begin{itemize}
\item
\tt{-DNDEBUG} will turn off debugging.  My recommendation is to leave debugging
on, since it doesn't slow things down too much.  Then when you decide to 
export a version of the code that needs to be as fast as possible, recompile
with this flag.  That said, the internal TMV code should be pretty well debugged
when you get it, so for compiling the library, you could go ahead and use this flag.

I keep two versions of the library for my own use - one with and one without 
this flag.  Then when I encounter a problem I can try linking to the debugging 
version to make sure it's not a bug in the TMV code.
\item
\tt{-DNWARN} will turn off warnings.  The code will output a few warnings to 
std::cout for a few things which aren't really errors, but are likely to be a
mistake in coding.  These can be turned off with this flag.
They are automatically turned off when -DNDEBUG is specified. 
\item
\tt{-DNOBLAS} will not call any external BLAS or LAPACK routines
\item
\tt{-DNOLAP} will not call any external LAPACK routines
\item
\tt{-DATLAS} will set up the BLAS calls as defined by ATLAS.  And (if -DNOLAP
is not specified), it will also call the several LAPACK routines
provided by ATLAS.
\item
\tt{-DCLAPACK} will set up the LAPACK calls for the CLAPACK distribution.  I 
find this version easier to get installed than the Fortran LAPACK distribution,
so I would recommend using this if you don't already have a version of LAPCK
installed somewhere on you system.  

Defining both ATLAS and CLAPACK will use
the CLAPACK version for all LAPACK routines, including the one also provided by 
ATLAS.  That is, ATLAS will only be used for its BLAS routines.  If you want
the ATLAS versions of its few LAPACK routines instead, the ATLAS 
installation instructions describe a way to get them into the CLAPACK library.
Also, you should make sure the INCLUDE specification lists the CLAPACK
version of clapack.h before the ATLAS version.
\item
\tt{-DMKL} will call all the external BLAS and LAPACK routines as defined by the
Intel Match Kernel Library.  You may want to check the file \tt{TMV\_Blas.h} to 
make sure the \tt{\#include} statement is correct for your installation of MKL.
\item
\tt{-DACML} will call all the external BLAS and LAPACK routines as defined by the 
AMD Core Math Library.  You may want to check the file \tt{TMV\_Blas.h} to make
sure the \tt{\#include} statement is correct for your installation of ACML.
\item
\tt{-DNO\_INST\_FLOAT} will not instantiate any \tt{<float>} classes or routines.
\item
\tt{-DNO\_INST\_DOUBLE} will not instantiate any \tt{<double>} classes or routines.
\item
\tt{-DNO\_INST\_COMPLEX} will not instantiate any complex classes or routines.
\item
\tt{-DINST\_LONGDOUBLE} will instantiate \tt{<long double>} classes and routines.
\item
\tt{-DINST\_INT} will instantiate \tt{<int>} classes and routines.
\item
\tt{-DNOSTL} uses some workarounds for segments of code which use the STL
library that didn't work for one of my compilers.  I'm pretty sure it is
because the compiler isn't completely installed correctly, so I don't think you
should really ever need to use this flag.  

But in any case, it 
will use a median-of-three quicksort algorithm for sorting
rather than the standard library's sort.  And it manually reads strings 
using character reads, rather than using just the $>>$ operator.  
\item
\tt{-DXTEST} will do extra testing in the test suite.  (It doesn't change anything
about the library.)  I always do my pre-release tests with this turned on, but the
executable gets quite large, as do many of the \tt{TMV\_Test*.o} files.  So I 
turn it off for the release version.  \tt{XTEST} mostly does additional checks
of the algorithms for different specific matrix pairs in the various
arithmetic operations to test possible failure modes.  If the shorter test suite
(i.e. without \tt{XTEST}) works ok for you, you should be fine.
\item
\tt{-DTMVFLDEBUG} will do extra (slow) debugging.  Specifically, 
it checks to make sure
element access is always within the range of the memory that was allocated.  
This shouldn't ever be necessary for a released version, since I should have already
discovered errors like this before the official release.  
But if you get strange segmentation
faults, you could give it a try to help me fix the problem.
\item
\tt{-DXDEBUG} will do different extra (even slower) debugging.  This one checks for
incorrect results from the various algorithms by doing things the simple slow way
and comparing the results to the fast blocked or recursive or in-place version
to make sure the answer isn't (significantly) different.  
I use this one a lot when debugging
new algorithms, usually on a file-by-file basis.  Again, you shouldn't need this for
an official release version.  But if you do get wrong answers for something, you could
use this to help me fix the problem.
\item
\tt{-DTMV\_BLOCKSIZE=NN} will change the block size used by some routines.  
Unless you have a strange CPU, I doubt this will speed things up much.  
But feel free to play around with it if you want.
\end{itemize}

\item{(advanced usage) Edit Inst and Blas files}

Note that, by default, the library will include instantiations of all classes and 
functions
which use \tt{<double>} or \tt{<float>} (including their complex versions).  
There are a few flags, such as 
\tt{-DNO\_INST\_FLOAT} and \tt{-DINST\_LONGDOUBLE},
which can change this as described above.  
But if you want to compile routines for some other class, say some user-defined
\tt{Quad} class, then you will need to modify the file \tt{TMV\_Inst.h}.  
You simply need to add the lines:
\begin{verbatim}
#define T Quad
#include InstFile
#undef T
\end{verbatim}
to the end of the file before compiling.  (Obviously, replace \tt{Quad} with whatever
type you want to be instantiated in the library.)

Also, the file \tt{TMV\_Blas.h} sets up all the BLAS and LAPACK calling structures,
as well as the necessary \tt{\#include} statements.  
So if the BLAS or LAPACK options aren't working 
for your system, you may need to edit these files as well.
This is especially true if your BLAS or LAPACK versions are not one of 
ATLAS, CLAPACK, MKL, or ACML.  The comments at the beginning of \tt{TMV\_Blas.h}
gives instructions on how to set up the file for other installations.

\item
Type: \\
\tt{make libs}

This will make the TMV libraries, \tt{libtmv.a} and \tt{libtmv\_symband.a},
which will be located in the directory lib.

Warning: this step may take a long time to finish.  So plan on letting this run
for an hour or so (depending on your machine and compiler of course).
Don't worry - once the library is compiled, compiling your code which uses the 
library is pretty quick.

I have tested the code using g++ versions 3.4.4 and 4.0.0, icc versions
8.0 and 9.0, and pgCC versions 5.2 and 6.1.  It should work with any ansi-compliant
compiler, but no guarantees if you use one other than these. 

Note: g++ version 3.3 has problems with some of the language features that
I use, so if you have this version or an older g++, you should probably
upgrade to a more recent version.

\item
(optional) Next type:\\
\tt{make tests}

This will make three executables called \tt{tmvtest1}, 
\tt{tmvtest2} and \tt{tmvtest3} in the bin directory.
This step also can take quite a while to finish, since the test code has tests
for lots of different combinations of matrix types in all the various 
arithmetic operators.

Then you should run the three test suites.
They should output a bunch of lines reading \tt{[{\em Something}] passed all tests}.
If one of them ends in a line that starts with \tt{Error}, 
then please send and email to \myemail\ about the problem.

You may also want to make the example programs by typing:\\
\tt{make examples}

This will make five executables called \tt{vector}, \tt{matrix}, \tt{division},
\tt{bandmatrix}, and \tt{symmatrix}.
These programs, along with their corresponding source code in the \tt{examples}
directory, give concrete examples of some of the common things you might want
to do with the TMV library.  They don't really try to be comprehensive,
but they do give a pretty good overview of the main features,
so looking at them may be a useful way to get started.

\item
Compile you program

Each .cpp file that uses TMV will need to have \\
\tt{\#include "TMV.h"}\\
at the top.
If you are using some of the special matrices, then you need to include their 
.h files as well.

To compile, you will need to use the compile flag
\tt{-I[tmvdir]/include} when making the object file to tell the 
compiler where the TMV header files are.

For the linking step, you need to compile with the flags
\tt{-L[tmvdir]/lib -ltmv -lm}.
If you are using any of \tt{SymMatrix}, \tt{HermMatrix}, \tt{BandMatrix}, 
\tt{SymBandMatrix}, or \tt{HermBandMatrix} in your code, then you will 
need to link with the flags \tt{-L[tmvdir]/lib -ltmv\_symband -ltmv -lm}.

If you are using BLAS and/or LAPACK calls from the TMV code,
then you will also need to 
link with their libraries.  
For example, for my version of Intel's MKL version of the BLAS and LAPACK
routines, I use 
\tt{-lmkl\_lapack -lmkl\_ia32 -lguide -lpthread}.  
For ATLAS,
I use \tt{-llapack -lcblas -latlas}.  For your specific installation, you may need 
the same thing, or something slightly different.

\end{enumerate}

Known compiler issues:
\begin{itemize}
\item
Apple's perversion of gcc does not work for compilation of the TMV library.  If you
want to compile this on a Mac, you should download the real gcc instead.  I 
recommend using Fink (http://fink.sourceforge.net/).

\item
The library is pretty big, so it can take quite a lot of memory to compile. 
For most compilers, it seems that a minimum of 256K to 512K is required.
(For compiling the test suite with the \tt{-DXTEST} flag, more than 
1GB of memory is recommended.)

\item
The Apple linker can't handle the size of the test suite's executable which compiled 
with \tt{-DXTEST}, even with plenty of memory.  So on this platform, you're 
pretty much stuck with just the regular test suites.  Even without \tt{-DXTEST}
you might need to compile some of the smaller test suites.  All of the tests
in \tt{tmvtest1} can be found in \tt{tmvtest1a}, \tt{tmvtest1b} and \tt{tmvtest1c}.
Likewise \tt{tmvtest2} has \tt{a}, \tt{b} and \tt{c} versions.  These are 
compiled by typing \tt{make test1a}, \tt{make test1b}, etc.  So if you want
to run the tests on a machine that can't link the full programs, these
smaller versions can help.

\item
The Portland Group C Compiler (or at least my installation of it) seems to 
do something
strange with the standard library's sort function. 

On one computer with pgCC version 5.2, the linker complains 
that it can't find the code for sort.  But sort should be completely inlined,
so the code should already be in \tt{TMV\_Vector.o}, where I use it.
Similarly, it has trouble linking the string read commands needed for reading
in \tt{SymBandMatrix}es.  Again, this should be included in the object file where 
I use it, \tt{TMV\_SymBandMatrix.o}.

On a different computer that has pgCC version 6.1, the sort command compiled
fine, but then didn't work correctly when I tested it.

The problems are probably due to something not being installed correctly on these
computers, or maybe I am just not including the correct linkages or something.  
But rather than trying to get the sysadmins of the
computers to find and fix the problem, 
I just added an option to compile with a simple median-of-three
quicksort algorithm, rather than the STL \tt{sort} command, and to read the strings
in character by character.  You can use this option by compiling with the 
flag \tt{-DNOSTL}.  

\item
If you get warnings (with compilations that don't use one of the 
\tt{-DNDEBUG} or \tt{-DNWARN} flags)
that look something like:
\begin{verbatim}
dgeqrf requested more workspace than provided 
for matrix with m,n = 1000,20 
Given: 2560, requested 5120
\end{verbatim}
then this means that your LAPACK uses a larger block size than 64.  This is 
a variable which is defined in \tt{TMV\_Blas.h} as \tt{LAP\_BLOCKSIZE}.
You should increase this to something larger.  In the above example, since
5120 (the ``requested'' size) is twice 2560 (the ``Given'' size), 
\tt{LAP\_BLOCKSIZE} should probably be increased by a factor of 2, so 128.
For your particular warning message, use the corresponding ratio (or larger)
for how much to increase \tt{LAP\_BLOCKSIZE}.
\end{itemize}

\newpage
\section{Known Bugs and Deficiencies (aka To Do List)}

If you find something to add to this list, or if you want me to bump something
to the top of the list, let me know.  Not that the list is currently in any kind of 
priority order, but, you know what I mean.  Email me at \myemail.

\begin{enumerate}
\item
\textbf{Symmetric arithmetic}

When writing complicated equations involving symmetric or hermitian matrices, 
you may find that an equation that seems perfectly ok does not compile.
The reason for this problem is explained in \S\ref{symarith} in some detail, 
so you should read about it there.  But basically, the workaround is usually
to break your equation up into smaller steps which do not require the code to 
explicitly instantiate any matrices.  For example: (this is the example from \S\ref{symarith})
\begin{verbatim}
s3 += x*s1 + s2;
\end{verbatim}
will not compile if \tt{s1}, \tt{s2}, and \tt{s3} are all symmetric, even though it is 
valid, mathematically.  Rewriting this as:
\begin{verbatim}
s3 += x*s1;
s3 += s2;
\end{verbatim}
will compile and work correctly.

\item
\textbf{Eigenvalues and eigenvectors}

I still need to figure out what interface I want for this to be as 
intuitive as possible.  Probably have a contained object which
contains the eigenvalue stuff similar to the way I keep
the decompositions for dividing.  

The code does calculate eigenvalues and eigenvectors
of hermitian matrices already, but the interface isn't intuitive.
The SVD of a hermitian matrix \underline{is} the 
eigenvector decomposition.  Well, almost.  Technically,
the singular values are defined to be positive, whereas the 
eigenvalues may be negative.  But when I calculate the SVD
for a \tt{HermMatrix},
I don't bother making them positive, since I figured people might
want to access the actual eigenvalues.  The columns of 
$U$ are the corresponding eigenvectors.

The code doesn't currently have any mechanism for calculating
eigenvalues for non-hermitian matrices.

\item
\textbf{More special matrix varieties}

Block-diagonal, generic sparse (probably both
row-based and column-based), block sparse, symmetric and hermitian block
diagonal, small symmetric and hermitian...

Maybe skew-symmetric and skew-hermitian.  Are these worth adding?

\item
\textbf{Packed storage}

Triangle and symmetric matrices. can be stored in (approximately) half the 
memory as a full $N \times N$ matrix using what is known as packed storage.  
There are BLAS routines for dealing with
these packed storage matrices, but I don't yet have the ability to 
create/use such matrices.

\item
\textbf{Small Specializations}

The SmallMatrix and SmallVector algorithms are all done inline which will
probably be close to optimally fast for small matrices, but I should look into
whether further algorithmic specialization for particular N,M values
would be worthwhile.

It may also be worth specializing a few of the divisions inline.  In particular,
I suspect division by a $2 \times 2$ matrix would be worth doing.

\item
\textbf{SVD algorithm}

My calculation of the 
SVD of a bidiagonal matrix uses the relatively slow QR algorithm,
rather than the divide-and-conquer or ``Relatively Robust Representation''
algorithms, which LAPACK uses.

I think these are the only routines for which the LAPACK version is still much
faster than my own.  The difference can be more than a factor of 10 for large 
matrices.

\item
\textbf{Row-major Bunch-Kauffman}

The Bunch-Kauffman decomposition for row-major symmetric/hermitian
matrices is currently $L D L^\dagger$, rather than $L^\dagger D L$.  
The latter should be somewhat (30\%?) faster.  The current $L D L^\dagger$
algorithm is the faster algorithm for column-major matrices.  
(These comments hold when the storage of the symmetric matrix is in the 
lower triangle - it is the opposite for storage in the upper triangle.)

\item
\textbf{Faster MultMV, MultMM}

I think I can mimic some of the ATLAS-style structure for the 
matrix-vector product and matrix-matrix product calculations 
to speed up the code for non-BLAS calculations.  

\item
\textbf{Conditions}

Currently, the SVD is the only decomposition which calculates the condition
of a matrix (specifically, the 2-condition).  
LAPACK has routines to calculate the 1- and infinity-condition
from an LU decomposition (and others).  I should add a similar capability.

\item
\textbf{Division error estimates}

LAPACK provides these.  It would be nice to add something along the same lines.

\item
\textbf{Equilibrate matrices}

LAPACK can equilibrate matrices before division.  Again, I should include this
feature too.  Probably as an option (since most matrices don't really need it)
as something like \tt{m.Equilibrate()} before calling a division routine.

\item
\textbf{Config program}

Currently, the user must modify the \tt{Makefile} and possibly 
\tt{TMV\_Blas.h}.  Most packages use some kind of configuration program
that asks the installer some questions and automatically sets everything
up correctly.  I've never written one of these, so this is a bit 
daunting for me, but it would make the 
installation a lot easier for users of the library.

\end{enumerate}

\newpage
\section{History}
\label{history}

\begin{description}
\item[Version 0.1]
The first matrix/vector library I wrote.  It wasn't very good, really. 
It had a lot of the functionality I needed, like mixing complex/real, 
SV decomposition, LU decomposition, etc.  But it wasn't at all fast for 
large matrices.  It didn't call BLAS or LAPACK, nor were the internal routines
very well optimized.  Also, while it had vector views for rows and columns, 
it didn't have matrix views for things like transpose.  Nor did it have any
delayed arithmetic evaluation.  And there were no special matrices.

I didn't actually name this one 0.1 until I had what I called version 0.3.

\item[Version 0.2]
This was also not named version 0.2 until after the fact.  It now had most of 
the current interface for regular Matrix and Vector operations.  I added
Upper/Lower TriMatrix and DiagMatrix.  I also had matrix views and
matrix composites to delay arithmetic evaluation.  The main problem was that
it was still slow.  I hadn't included any BLAS calls yet.  And while the 
internal routines at least used algorithms that used unit strides whenever 
possible, they didn't do any blocking or recursion, which are key for large
matrices.

\item[Version 0.3]
Finally, I actually named this one 0.3 at the time.  (Actually, at first and also 
for 0.2, I was calling it JMV rather than TMV for ``Jarvis's Matrix/Vector library''.
Then I decided that it was too egotistical, so I switched the first letter to 
T for ``Template''.)  The big addition here was BLAS and LAPACK calls, which 
helped me to realize how slow my internal code really was (although I hadn't
updated them to block or recursive algorithms yet).

I also added BandMatrix.

\item[Version 0.4]
The new version number here was because I needed some added functionality for a
project I was working on.  It retrospect, it really only deserves a 0.01 increment,
since the changes weren't that big.  But, oh well.

\begin{itemize}
\item
Added QR\_Downdate.  (This was the main new functionality I needed.)
\item
Improved the numerical accuracy of the QRP decomposition.
\item
Added the possibility of not storing U,V for the SVD.
\item
Greatly improved the test suite, and correspondingly found and corrected a few bugs.
\item
Added some arithmetic functionality that had been missing (like \tt{m += L*U}).
\end{itemize}

\item[Version 0.5]
The new symmetric matrix classes precipitated a major version number update.
I also sped up a lot of the algorithms:
\begin{itemize}
\item
Added SymMatrix, HermMatrix, and all associated functionality.
\item
Finally added blocked versions of most of the algorithms, so the non-LAPACK
code runs a lot faster.
\item
Allowed for loose QRP decomposition.
\item
Added the possibility of division in place, rather than copying to new storage.
\end{itemize}

\item[Version 0.51]
Some minor improvements:
\begin{itemize}
\item
Sped up some functions like matrix addition and assignment by adding the 
LinearView method, which can view a matrix a one long vector.
\item
Added QR\_Update.
\item
Blocked some more algorithms like TriMatrix multiplication/division, so non-BLAS
code runs significantly faster (but still slower than BLAS).
\end{itemize}

\item[Version 0.52]
The first ``public'' release!  And correspondingly, the first with documentation
and a web site.  A few other people had used previous versions, but since the
only documentation was my comments in the .h files, it wasn't all that user-friendly.
\begin{itemize}
\item
Added SaveDiv() and related methods for division control.  It used to be that the 
default behavior was to reuse the matrix decomposition once it was made.  
This meant that naive code could get the wrong answer for routines which
divide, change the matrix, and then divide again.  (ReSetDiv() had been available
to redo the decomposition, but you had to remember to use it.)  I decided
(on the advice of Gary Bernstein) that it was better to have naive code 
be correct and slow, rather than fast and wrong.  So now the programmer needs 
to explicitly say to save the decomposition.
\item
Extended the test suite a bit and found a couple more (fairly obscure) bugs.
Mostly with non-square or very fat band matrices.
\item
Added in-place versions of the algorithms for $S = L^\dagger L$ and $S = L L^\dagger$.
Basically, this is Cholesky in reverse.
\end{itemize}

\item[Version 0.53]
By popular demand (well, a request by Fritz Stabenau, at least):
\begin{itemize}
\item
Added the Fortran-style indexing.
\end{itemize}

\item[Version 0.54]
Inspired by my to-do list, which I wrote for Version 0.52, I tackled a few of the 
items on the list and addressed some issues that people had been having
with compiling:
\begin{itemize}
\item
Changed from a rudimentary exception scheme (with just one class - 
\tt{tmv\_exception}) to the current more full-featured exception hierarchy.
Also added \tt{auto\_ptr} usage instead of bald pointers to make the 
exception throws memory-leak safe.
\item
Sped up SymLUDiv and SymSVDiv inverses.
\item
Added the possibility of compiling a small version of the library and test suite.
\item
Fixed some things which gave warnings by pgCC and by icc with the -Wall flag
(a far more extensive list of warnings than gcc's -Wall).  There are still a 
few things that icc -Wall warns about that I'm not going to try to fix.  
Namely remarks numbered 279, 383, 810, and 981.  The code passes all the rest
of their warnings and remarks.
\item
Consolidated SymLUDiv and HermLUDiv classes into just SymLUDiv, which now checks
whether the matrix is hermitian automatically.  Unlike with SymSVDiv and
HermSVDiv, there is no good reason to have these be separate classes.
\item
Reduced the number of operations that make temporary matrices when multiple
objects use the same storage.  For example, $M = L U$ can be done in place now, where 
$L$ and $U$ are the lower and upper triangle portions of $M$.
\item
Fixed some syntax problems with symmetric matrix arithmetic.
\item
Specialized Tridiagonal $\times$ Matrix calculation.
\item
Added ElementProd and AddElementProd functions for matrices (rather than just
for vectors).
\item
Added I/O tests to the test suite, and found/fixed a couple of bugs in these 
routines.
\item
Added CLAPACK and ACML as known versions of LAPACK (the LAPACK calls
had previously been specific to the Intel MKL version).  Also made the file
\tt{TMV\_Blas.h} much better organized, so it is easier to tailor the code to 
a different LAPACK distribution with different calling conventions.

\end{itemize}

\item[Version 0.60]
This revision merits a first-decimal-place increment, since I added a 
few big features.  I also registered it with SourceForge, 
which is a pretty big milestone as well.
\begin{itemize}
\item
Added \tt{SmallVector} and \tt{SmallMatrix} with all accompanying algorithms.
\item 
Added \tt{SymBandMatrix} and \tt{HermBandMatrix} with all accompanying algorithms.
\item
Made arithmetic between any two special matrices compatible, so long as 
the operation is allowed given their respective shapes.  
e.g. \tt{U += B} and \tt{U *= B} are allowed if \tt{B} is upper-banded.  
There are many other examples which are now legal statements.
\item
Changed \tt{QR\_Downdate()} to throw an exception rather than return false
when it fails.
\item
Added the GPL License to the end of this documentation, and a copyright
and GPL snippet into each file.
\item
Changed the -D compiler options for changing which types get instantiated.
The default is still to instantiate double and float.  Now to turn these
off, use -DNO\_INST\_DOUBLE and -DNO\_INST\_FLOAT respectively.  To
add int or long double instantiations, use -DINST\_INT and 
-DINST\_LONGDOUBLE respectively.
\item
Split up the library into \tt{libtmv.a} and \tt{libtmv\_symband.a}.  The
latter includes everything for the \tt{SymMatrix}, \tt{HermMatrix}, 
\tt{BandMatrix}, \tt{SymBandMatrix}, and \tt{HermBandMatrix} classes.
The former includes everything else including \tt{DiagMatrix} and
\tt{Upper/LowerTriMatrix}.  Also, I got rid of \tt{libsmalltmv.a}, which
was basically the same as the new trimmed down \tt{libtmv.a}.
\item 
Found and fixed various minor bugs along the way.

\end{itemize}

\end{description}

\newpage
\section{License}
\label{license}

This software is licensed under the GNU General Public License.  This license is
sometimes referred to as ``copyleft''.  This gist of the license (made more specific
in the following text) is that you are allowed to copy, modify, and use TMV code
with basically no restrictions as far as your own use is concerned, 
but with a few restrictions regarding further distribution.

The biggest restriction is that if you distribute code which
uses the TMV code, then you must distribute it with a copyleft license as well.
In other words, you may not use the TMV library with commercial, proprietary 
software. 

The other main restriction on such a distribution is that you must include this 
document,
including the below license, along with your version of the code.  If you modified
any portions of the library, you must indicate your modifications in the accompanying
documentation files.

More information may be obtained from http://www.gnu.org/copyleft/.

\begin{center}
{\parindent 0in

Version 2, June 1991

Copyright \copyright\ 1989, 1991 Free Software Foundation, Inc.

\bigskip

51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA

\bigskip

Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
}
\end{center}

\renewcommand{\abstractname}{Preamble}
\begin{abstract}
The licenses for most software are designed to take away your freedom to
share and change it.  By contrast, the GNU General Public License is
intended to guarantee your freedom to share and change free software---to
make sure the software is free for all its users.  This General Public
License applies to most of the Free Software Foundation's software and to
any other program whose authors commit to using it.  (Some other Free
Software Foundation software is covered by the GNU Library General Public
License instead.)  You can apply it to your programs, too.

When we speak of free software, we are referring to freedom, not price.
Our General Public Licenses are designed to make sure that you have the
freedom to distribute copies of free software (and charge for this service
if you wish), that you receive source code or can get it if you want it,
that you can change the software or use pieces of it in new free programs;
and that you know you can do these things.

To protect your rights, we need to make restrictions that forbid anyone to
deny you these rights or to ask you to surrender the rights.  These
restrictions translate to certain responsibilities for you if you
distribute copies of the software, or if you modify it.

For example, if you distribute copies of such a program, whether gratis or
for a fee, you must give the recipients all the rights that you have.  You
must make sure that they, too, receive or can get the source code.  And
you must show them these terms so they know their rights.

We protect your rights with two steps: (1) copyright the software, and (2)
offer you this license which gives you legal permission to copy,
distribute and/or modify the software.

Also, for each author's protection and ours, we want to make certain that
everyone understands that there is no warranty for this free software.  If
the software is modified by someone else and passed on, we want its
recipients to know that what they have is not the original, so that any
problems introduced by others will not reflect on the original authors'
reputations.

Finally, any free program is threatened constantly by software patents.
We wish to avoid the danger that redistributors of a free program will
individually obtain patent licenses, in effect making the program
proprietary.  To prevent this, we have made it clear that any patent must
be licensed for everyone's free use or not licensed at all.

The precise terms and conditions for copying, distribution and
modification follow.
\end{abstract}

\begin{center}
{\Large \sc GNU General Public License
\\\vspace{3mm}Terms and Conditions For Copying, Distribution and Modification}
\end{center}


\begin{enumerate}

\addtocounter{enumi}{-1}

\item 

This License applies to any program or other work which contains a notice
placed by the copyright holder saying it may be distributed under the
terms of this General Public License.  The ``Program'', below, refers to
any such program or work, and a ``work based on the Program'' means either
the Program or any derivative work under copyright law: that is to say, a
work containing the Program or a portion of it, either verbatim or with
modifications and/or translated into another language.  (Hereinafter,
translation is included without limitation in the term ``modification''.)
Each licensee is addressed as ``you''.

Activities other than copying, distribution and modification are not
covered by this License; they are outside its scope.  The act of
running the Program is not restricted, and the output from the Program
is covered only if its contents constitute a work based on the
Program (independent of having been made by running the Program).
Whether that is true depends on what the Program does.

\item You may copy and distribute verbatim copies of the Program's source
  code as you receive it, in any medium, provided that you conspicuously
  and appropriately publish on each copy an appropriate copyright notice
  and disclaimer of warranty; keep intact all the notices that refer to
  this License and to the absence of any warranty; and give any other
  recipients of the Program a copy of this License along with the Program.

You may charge a fee for the physical act of transferring a copy, and you
may at your option offer warranty protection in exchange for a fee.

\item

You may modify your copy or copies of the Program or any portion
of it, thus forming a work based on the Program, and copy and
distribute such modifications or work under the terms of Section 1
above, provided that you also meet all of these conditions:

\begin{enumerate}

\item 

You must cause the modified files to carry prominent notices stating that
you changed the files and the date of any change.

\item

You must cause any work that you distribute or publish, that in
whole or in part contains or is derived from the Program or any
part thereof, to be licensed as a whole at no charge to all third
parties under the terms of this License.

\item
If the modified program normally reads commands interactively
when run, you must cause it, when started running for such
interactive use in the most ordinary way, to print or display an
announcement including an appropriate copyright notice and a
notice that there is no warranty (or else, saying that you provide
a warranty) and that users may redistribute the program under
these conditions, and telling the user how to view a copy of this
License.  (Exception: if the Program itself is interactive but
does not normally print such an announcement, your work based on
the Program is not required to print an announcement.)

\end{enumerate}


These requirements apply to the modified work as a whole.  If
identifiable sections of that work are not derived from the Program,
and can be reasonably considered independent and separate works in
themselves, then this License, and its terms, do not apply to those
sections when you distribute them as separate works.  But when you
distribute the same sections as part of a whole which is a work based
on the Program, the distribution of the whole must be on the terms of
this License, whose permissions for other licensees extend to the
entire whole, and thus to each and every part regardless of who wrote it.

Thus, it is not the intent of this section to claim rights or contest
your rights to work written entirely by you; rather, the intent is to
exercise the right to control the distribution of derivative or
collective works based on the Program.

In addition, mere aggregation of another work not based on the Program
with the Program (or with a work based on the Program) on a volume of
a storage or distribution medium does not bring the other work under
the scope of this License.

\item
You may copy and distribute the Program (or a work based on it,
under Section 2) in object code or executable form under the terms of
Sections 1 and 2 above provided that you also do one of the following:

\begin{enumerate}

\item

Accompany it with the complete corresponding machine-readable
source code, which must be distributed under the terms of Sections
1 and 2 above on a medium customarily used for software interchange; or,

\item

Accompany it with a written offer, valid for at least three
years, to give any third party, for a charge no more than your
cost of physically performing source distribution, a complete
machine-readable copy of the corresponding source code, to be
distributed under the terms of Sections 1 and 2 above on a medium
customarily used for software interchange; or,

\item

Accompany it with the information you received as to the offer
to distribute corresponding source code.  (This alternative is
allowed only for noncommercial distribution and only if you
received the program in object code or executable form with such
an offer, in accord with Subsection b above.)

\end{enumerate}


The source code for a work means the preferred form of the work for
making modifications to it.  For an executable work, complete source
code means all the source code for all modules it contains, plus any
associated interface definition files, plus the scripts used to
control compilation and installation of the executable.  However, as a
special exception, the source code distributed need not include
anything that is normally distributed (in either source or binary
form) with the major components (compiler, kernel, and so on) of the
operating system on which the executable runs, unless that component
itself accompanies the executable.

If distribution of executable or object code is made by offering
access to copy from a designated place, then offering equivalent
access to copy the source code from the same place counts as
distribution of the source code, even though third parties are not
compelled to copy the source along with the object code.

\item
You may not copy, modify, sublicense, or distribute the Program
except as expressly provided under this License.  Any attempt
otherwise to copy, modify, sublicense or distribute the Program is
void, and will automatically terminate your rights under this License.
However, parties who have received copies, or rights, from you under
this License will not have their licenses terminated so long as such
parties remain in full compliance.

\item
You are not required to accept this License, since you have not
signed it.  However, nothing else grants you permission to modify or
distribute the Program or its derivative works.  These actions are
prohibited by law if you do not accept this License.  Therefore, by
modifying or distributing the Program (or any work based on the
Program), you indicate your acceptance of this License to do so, and
all its terms and conditions for copying, distributing or modifying
the Program or works based on it.

\item
Each time you redistribute the Program (or any work based on the
Program), the recipient automatically receives a license from the
original licensor to copy, distribute or modify the Program subject to
these terms and conditions.  You may not impose any further
restrictions on the recipients' exercise of the rights granted herein.
You are not responsible for enforcing compliance by third parties to
this License.

\item
If, as a consequence of a court judgment or allegation of patent
infringement or for any other reason (not limited to patent issues),
conditions are imposed on you (whether by court order, agreement or
otherwise) that contradict the conditions of this License, they do not
excuse you from the conditions of this License.  If you cannot
distribute so as to satisfy simultaneously your obligations under this
License and any other pertinent obligations, then as a consequence you
may not distribute the Program at all.  For example, if a patent
license would not permit royalty-free redistribution of the Program by
all those who receive copies directly or indirectly through you, then
the only way you could satisfy both it and this License would be to
refrain entirely from distribution of the Program.

If any portion of this section is held invalid or unenforceable under
any particular circumstance, the balance of the section is intended to
apply and the section as a whole is intended to apply in other
circumstances.

It is not the purpose of this section to induce you to infringe any
patents or other property right claims or to contest validity of any
such claims; this section has the sole purpose of protecting the
integrity of the free software distribution system, which is
implemented by public license practices.  Many people have made
generous contributions to the wide range of software distributed
through that system in reliance on consistent application of that
system; it is up to the author/donor to decide if he or she is willing
to distribute software through any other system and a licensee cannot
impose that choice.

This section is intended to make thoroughly clear what is believed to
be a consequence of the rest of this License.

\item
If the distribution and/or use of the Program is restricted in
certain countries either by patents or by copyrighted interfaces, the
original copyright holder who places the Program under this License
may add an explicit geographical distribution limitation excluding
those countries, so that distribution is permitted only in or among
countries not thus excluded.  In such case, this License incorporates
the limitation as if written in the body of this License.

\item
The Free Software Foundation may publish revised and/or new versions
of the General Public License from time to time.  Such new versions will
be similar in spirit to the present version, but may differ in detail to
address new problems or concerns.

Each version is given a distinguishing version number.  If the Program
specifies a version number of this License which applies to it and ``any
later version'', you have the option of following the terms and conditions
either of that version or of any later version published by the Free
Software Foundation.  If the Program does not specify a version number of
this License, you may choose any version ever published by the Free Software
Foundation.

\item
If you wish to incorporate parts of the Program into other free
programs whose distribution conditions are different, write to the author
to ask for permission.  For software which is copyrighted by the Free
Software Foundation, write to the Free Software Foundation; we sometimes
make exceptions for this.  Our decision will be guided by the two goals
of preserving the free status of all derivatives of our free software and
of promoting the sharing and reuse of software generally.

\begin{center}
{\Large\sc
No Warranty
}
\end{center}

\item
{\sc Because the program is licensed free of charge, there is no warranty
for the program, to the extent permitted by applicable law.  Except when
otherwise stated in writing the copyright holders and/or other parties
provide the program ``as is'' without warranty of any kind, either expressed
or implied, including, but not limited to, the implied warranties of
merchantability and fitness for a particular purpose.  The entire risk as
to the quality and performance of the program is with you.  Should the
program prove defective, you assume the cost of all necessary servicing,
repair or correction.}

\item
{\sc In no event unless required by applicable law or agreed to in writing
will any copyright holder, or any other party who may modify and/or
redistribute the program as permitted above, be liable to you for damages,
including any general, special, incidental or consequential damages arising
out of the use or inability to use the program (including but not limited
to loss of data or data being rendered inaccurate or losses sustained by
you or third parties or a failure of the program to operate with any other
programs), even if such holder or other party has been advised of the
possibility of such damages.}

\end{enumerate}


\begin{center}
{\Large\sc End of Terms and Conditions}
\end{center}

This library is largely based upon work supported by the 
National Aeronautics and Space Administration (NASA) under Grant No 
NAG5-10924, a Keck Foundation grant, and the University of Pennsylvania. 
Any opinions, findings and conclusions or recommendations expressed in this 
material are those of the author(s) and do not necessarily reflect the views 
of NASA, the Keck Foundation, or the University of Pennsylvania.
This software is made available to you on an ``as is'' basis with no 
representations or warranties, express or implied, including but not 
limited to any warranty of performance, merchantability, fitness for a 
particular purpose, commercial utility, non-infringement or title.  
Neither the author nor the organizations providing the support under 
which the work was developed will be liable to you or any third party 
with respect to any claim arising from your further development of the 
software or any products related to or derived from the software, or for 
lost profits, business interruption, or indirect special or consequential 
damages of any kind. 

\end{document}
